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1.0  INTRODUCTION 

Rapid,  time-resolved  diagnostics  of  various  particle-laden  flows  are  of  much 
importance  in  specific  Air  Force  applications.  In  aerodynamic  simulation  facilities,  the 
presence  of  dust,  dirt,  erosion  products,  and  condensate  can  adversely  affect  both  the 
operation  of  the  test  facilities  and  the  interpretation  of  measurements  made  in  those 
facilities.  For  such  activities,  knowledge  of  the  particle  size  distribution  function  (PSDF) 
is  valuable.  Occasionally,  the  presence  of  natural  particulate  matter  in  a  test  facility  is 
desirable,  as  is  the  case  when  using  laser  velocimetry  to  determine  the  velocity  in 
unseeded  flow  fields.  Even  then,  a  knowledge  of  the  PSDF  is  required  to  characterize  the 
gas-particle  interaction  and,  especially,  to  indicate  whether  the  particles  being  observed 
can  be  expected  to  follow  the  flow  sufficiently  well  to  permit  accurate  measurements. 

Rocket  exhaust  diagnostics  and  spacecraft  contamination  studies  also  require  a 
known  PSDF  and  particle  number  density  within  local  flow  regions.  In  the  case  of  rocket 
exhausts,  it  is  desirable  to  know  the  total  number  density,  chemical  species  distribution, 
and  the  PSDF  of  the  particulate  effluent  (l)  to  understand  primary  and  secondary  smoke 
formation  and  to  validate  prediction  models  of  these  phenomena  and  (2)  to  characterize 
the  particulate  exhaust  effluents  of  both  liquid-  and  solid-propellant  rockets.  Such 
characterization  is  required  to  assess  the  contaminant  potential  of  such  effluents. 

Also,  the  diagnostic  techniques  discussed  in  this  report  provide  a  unique  means  of 
studying  one  of  the  most  important  and  least  understood  aspects  of  two-phase  flow  fields 
—  the  formation  of  liquid  droplets  and  small  solid  particles  from  the  primitive  molecular 
clusters  (dimers,  trimers,  tetramers,  and  so  forth)  which  can  sometimes  result  in  massive 
condensation  in  hypersonic  flow  expansions. 

Before  the  development  program  discussed  in  this  report  was  begun,  it  was  decided 
that  the  proposed  diagnostics  scheme  should  meet  the  following  requirements: 

1.  The  technique  should  not  be  restricted  to  the  requirement  that  a  single 
particle  be  present  in  the  focal  volume,  but,  rather,  the  method  should  be 
capable  of  yielding  accurate  results  when  a  polydisperse  size  distribution  is 
present. 

2.  Both  average  particle  sizes  and  PSDF’s  should  be  obtainable  from  the 
technique,  and  submicron  particles  should  be  measurable. 

3.  Time  resolution  should  be  possible  to  within  milliseconds. 
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4.  It  should  be  possible  to  evaluate  uncertainties  in  the  measurements  and  to 
estimate  the  quantitative  effect  of  measurement  imprecision  on  the 
particle  size  results. 

5.  Insofar  as  possible,  the  technique  should  require  only  standard  laser 
sources  and  detectors  (rather  than  requiring  dramatic  advances  in  the  state 
of  the  art  of  optical  and  electronic  instrumentation). 

a 

6.  The  technique  should  not  require  absolute  intensity  measurements. 

7.  Although  originally  intended  to  size  particles  of  a  single  chemical  species 
with  known  complex  indices  of  refraction,  the  technique  should  be 
capable  of  generalization  (1)  to  the  measurement  of  the  PSDF  of  a  single 
species  of  unknown  refractive  index  and  (2)  to  the  determination  of 
effective  complex  refractive  indices  for  a  particle  mixture  whose  chemical 
identities  are  not  known  a  priori.  In  addition  to  the  latter,  the  technique 
should  permit  the  deconvolution  of  both  particle  species  and  size 
distributions. 

To  satisfy  these  requirements,  it  was  decided  that  the  program  of  particle  diagnostics 
should  be  based  on  multiangular  or  multispectral  Mie  scattering  measurements  or, 
perhaps,  on  a  combination  of  both  types  of  measurements.  Although  only  the 
multiangular  measurements  and  analyses  have  been  used,  the  combination  of 
measurements  may  well  be  required  by  some  future  applications.  Mie  scattering 
measurements  have  a  distinct  advantage  over  measurements  based  on  more  approximate 
analysis  schemes,  since  the  Mie  theory  is  an  exact  solution  for  the  polarized,  single 
scattering  of  light  by  a  dielectric  sphere  of  arbitrary  size.  However,  the  optical  depth  for 
Mie  scattering  must  be  small  enough  to  ensure  that  multiple-scattering  events  contribute 
no  more  than  from  approximately  10  to  20  percent  to  the  scattering  signal  at  any  given 
detector  location.  Also,  departure  from  sphericity  affects  the  accuracy  of  interpretation 
of  the  polarized  measurements  in  the  backscattering  hemisphere.  In  addition,  the 
refractive  indices  of  the  measured  particles  must  be  known  or  be  iteratively  determinable 
within  the  sizing  analysis  procedure.  Various  procedures  exist  to  approximate 
nonspherical  scatterers,  and  an  assessment  of  their  use,  in  conjunction  with  the  sizing 
procedures  discussed  in  this  report,  is  underway.  However,  a  complete  assessment  of  the 
effects  of  nonspherical  scattering  is  beyond  the  scope  of  this  report;  the  results  discussed 
pertain  only  to  scattering  from  spheres.  Subsequent  sections  of  this  report  describe  and 
summarize  the  analytical  procedures  that  have  been  developed  or  extended  at  AEDC 
for  Mie  scattering  PSDF  determinations. 
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1.1  REVIEW  OF  AEDC  MIE  SCATTERING  PARTICLE 
DIAGNOSTICS  APPLICATIONS 

The  development  of  the  particle  diagnostic  procedures  occurred  in  several  stages,  the 
first  of  which  was  the  development  of  a  comprehensive  Mie  angular  scattering  computer 
code  (Ref.  1)  that  was  based  on  routines  reported  in  Ref.  2.  To  be  applied  to  particle 
diagnostics  problems  at  AEDC,  the  program  had  to  be  modified  in  several  ways; 
consequently,  a  family  of  codes,  collectively  known  as  MIESTEP,  was  generated.  These 
codes  compute,  plot,  and  store  on  tape  files  the  Mie  scattering  functions,  which  are 
integrated  over  specified  angles,  and  ancillary  scattering  functions  for  a  specified  range  of 
size  parameters.  These  analytical  procedures  and  a  primitive  version  of  the  size 
deconvolution  code  discussed  in  Appendix  A  were  used  to  analyze  Mie  scattering  data 
obtained  in  the  particulate-laden  flow  of  an  arc-driven  hypervelocity  wind  tunnel  at 
AEDC.  These  results,  reported  in  Ref.  3  and  Appendix  B,  established  the  average  particle 
size,  a  four-bin  histogram  approximation  to  the  PSDF,  and  the  particle  number  density  of 
these  submicron  particles.  Also,  the  results  obtained  with  the  measurements  and  analysis 
procedure  were  used  to  eliminate  water  and  iron  as  possible  dominant  scattering  species 
and  to  establish  limits  for  the  mass  fractions  of  the  dominant  scattered,  which  were 
copper  and  tungsten  produced  by  the  arc.  Furthermore,  these  results  were  consistent  with 
a  posttest  chemical  analysis  of  electrode  weight  loss  (See  Appendix  B.) 

After  this  successful  application  of  Mie  scattering  principles  to  particle  diagnostics  in 
an  AEDC  test  facility,  a  series  of  laboratory  measurements  of  the  scattering  from  a 
suspension  of  polystyrene  spheres  was  made  at  both  single-  and  multiple-scattering 
conditions.  The  purpose  of  this  work  [performed  under  the  sponsorship  of  the  Air  Force 
Rocket  Propulsion  Laboratory  (AFRPL)  and  reported  in  Ref.  4]  was  to  validate  a 
multiple-scattering  code  that  was  developed  by  an  AFRPL  contractor,  under  conditions 
of  simulated  scattering  from  solid-propellant  rocket  motor  (SRM)  plumes.  Figure  1,  taken 
from  Ref.  4,  shows  the  agreement  of  Mie  theory  with  the  single-scattering,  monodisperse 
scattering  data  obtained  in  Ref.  4.  Figure  2  shows  a  similar  comparison  to  results  that 
were  reported  in  Ref.  5  and  that  were  obtained  from  a  narrow,  log-normal  distribution  of 
spherical  aerosols. 

Also,  a  series  of  flow  visualization  measurements  was  made  for  the  MIT  Lincoln 
Laboratory  under  the  sponsorship  of  the  Air  Force  Space  and  Missile  Systems 
Organization  (SAMSO).  These  experiments  (reported  in  Ref.  6)  involved  using  an  image 
intensifier  to  permit  high-speed  photography  of  the  flow  field  surrounding  particles 
ejected  from  a  model  into  a  hypervelocity  flow  field.  Although  Mie  scattering  sizing 
measurements  were  not  made,  Mie  calculations  based  on  approximate  and  independent 
pretest  size  measurements  permitted  the  quantitative  interpretation  of  isointensity 
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contours  (obtained  with  computerized  densitometry  of  the  photographs  of  the 
image-intensifier  surface)  as  contours  of  constant  particle  number  densities. 

One  of  the  procedures  was  applied  to  sizing  particles  with  diameters  of  about  0.1  n 
in  conjunction  with  tests  of  a  small  liquid-bipropellant  thruster.  An  image  intensifier, 
located  at  a  90-deg  scattering  angle,  was  used  for  flow  visualization  of  a  laser-illuminated 
exhaust  plume.  Analysis  of  the  degree-of-polarization  of  the  scattering  signal  yielded  a 
quantitative  measure  of  the  order  of  magnitude  of  the  particle  size  of  the  exhaust.  These 
results  were  reported  in  Ref.  7. 

The  particle  diagnostics  procedures  discussed  in  this  report  are  also  currently  being 
applied  in  the  BATES  motor  test  program  at  the  AFRPL.  These  measurements  of  the 
exhaust  plumes  of  SRM’s  are  being  made  with  various  diagnostics  techniques,  one  of 
which  is  a  Mie  scattering  determination  of  the  PSDF  involving  the  experimental 
configuration  shown  in  Fig.  3.  Data  obtained  from  these  experiments  will  be  analyzed 
using  the  size  deconvolution  procedures  presented  in  this  report.  Computer  simulations 
indicate  that  six-bin  histograms  of  the  PSDF  having  acceptable  accuracy  can  be  obtained 
if  the  imprecision  in  the  scattering  ratios  (which  are  inputs  to  the  deconvolution 
procedure)  can  be  restricted  to  average  values  (for  all  the  scattering  channels)  no  larger 
than  approximately  ±3  percent.  A  larger  input  imprecision  of,  for  example,  ±5  percent 
will  reduce  the  attainable  size  resolution  to  no  more  than  Five  size  bins,  which  span  the 
range  of  diameters  from  less  than  0.1  to  approximately  1.6  n. 

1.2  SUMMARY  OF  AEDC  MIE  SCATTERING  PARTICLE 
DIAGNOSTIC  METHODS 

Before  detailing  the  parameters  that  are  determined  by  the  AEDC  particle  diagnostic 
procedure,  it  is  worthwhile  to  summarize  the  present  form  of  the  particle  diagnostics 
approach  and  options;  Fig.  4  is  a  flow  chart  of  these  procedures.  Figure  4  shows  that 
several  parallel  paths  provide  a  degree  of  redundancy  in  particle  diagnostics. 
Measurements  made  at  several  scattering  angles  of  both  the  degree  of  polarization  and  of 
a  set  of  functions,  called  ANGER  functions  (Angular/Extinction  Ratios),  provide  a  means 
of  estimating  a  weighted  average  particle  size,  which,  if  the  PSDF  is  not  determined,  is  to 
be  interpreted  as  a  surface  mean  diameter. 

The  lower  half  of  the  flow  chart  indicates  the  procedures  used  for  deconvolving  the 
PSDF,  Input  scattering  ratios  (DECRATS)  are  used  to  deconvolve  the  PSDF  in  a 
histogram  representation.  The  computer  program  that  accomplishes  this  deconvolution  is 
detailed  in  Appendix  A.  Outputs  from  the  deconvolution  histogram  can  be  used  as  initial 
guesses  to  perform  regression  procedures  currently  under  development.  These  procedures 
will  permit  reconstructing  the  normalized  PSDF  as  one  of  the  following  three 
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single-mode,  parametric  size  distributions:  (1)  a  three-parameter,  Deirmendjian 
modified-gamma  distribution  (Ref.  8);  (2)  a  two-parameter,  log-normal  distribution;  and 
(3)  a  two-parameter,  truncated  Gaussian  distribution  that  is  normalized  on  a  semi-infinite 
size  interval.  If  the  size  deconvolution  histogram  indicates  the  occurrence  of  either  a 
bimodal  or  a  higher-order,  multimodal  distribution,  then  these  direct  regression 
procedures  are  generally  inappropriate. 

Determining  the  average  diameter  or  the  PSDF  permits  estimating  the  total  particle 
number  density  by  convolving  the  computed  extinction  cross-section  and  transmission 
measurements.  This  method  for  measuring  the  particle  number  density  is  discussed  in 
Ref.  3. 


2.0  THEORY 


2.1  BACKGROUND 


The  Mie  scattering  theory  is  detailed  in  many  textbooks  on  optics  or 
electromagnetic  theory;  Refs.  9  through  1 1  are  especially  useful.  To  understand  the 
properties  of  the  scattering  ratios  used  in  the  AEDC  particle  diagnostics  program,  one 
must  first  consider  the  application  of  Mie  theory  to  the  interpretation  of  light-scattering 
measurements. 


First,  one  should  consider  the  angular  distribution  of  single  scattering  of  polarized 
radiation  that  is  incident  on  a  monodisperse  distribution  of  particles  of  diameter  D.  If  the 
incident  beam  power  is  $OJ  then  flux,  which  is  scattered  with  polarization  vector 
perpendicular  to  the  plane  of  scattering  into  a  given  solid  angle  subtended  by  the 
detection  optics  scheme,  is  written  as  follows: 


^Tj  +  A^/2  j*  0.  +  A0./2 


r  ifr,  +  2  V. 

n(x,  n  /  / 


dfi 


(x,  6,  0,J?)  sin0  id  d£' 
(1) 


It  is  assumed  that  the  incident  beam  is  along  the  z-axis  and  that  the  incident  polarization 
direction  is  the  x-axis  (see  Fig.  5).  The  parallel  polarization  component  is  described  by 


i(/i  +  A^./2  -  0.  +  A0./ 2  dp j 


a .  +■  tip./  z  aw, 

- (x*  6 ,  0  i?)  sin#  d0  Atf/  AV 

e.  -  A0./2 

(2) 


In  these  equations  the  solid  angle,  Aft,  subtended  by  the  detector  at  the  scattering 
volume  of  length  £  is  given  by 
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Ail  =  sin  8  A0  A^r 

One  more  equation  completes  the  set  of  equations  that  describes  the  single 
scattering  and  attenuation  by  spheres*.  The  transmitted  flux  is  given  by 


*Tr  <*-  >?)/*„ 


0  <r„xl  (x,  7})  df 


] 


T(X,  7,) 


where  aext  is  the  total  extinction  cross  section 

°ext  1>  -  ^cat  +  ffabs  (x-  7> 


(3) 


and  L  is  the  total  distance  over  which  scattering  and  absorption  occur.  The  total 
scattering  cross  section  is  the  integral  of  the  differential  cross  section  for  scattering 
unpolarized  light  over  the  entire  spherical  solid  angle,  4 ir  sr;  i.e. , 

/■iff 

—  (x,  6,  ijj,  jj)  d£ 


where 


dffy  dOj  dff2 

dQ  _  dH  +  dfi 

The  absorption  cross  section  and  the  total  scattering  cross  section  are  the  products  of  the 
particle’s  geometric  cross  section,  7r(D/2)2,  and  the  dimensionless  scattering  and 
absorption  efficiency  factors,  Qscat  and  Qabs,  respectively.  These  two  analytical  functions 
are  obtained  directly  (as  functions  of  particle  size  and  complex  refractive  index)  from  the 
Mie  scattering  series  solution,  without  the  need  for  numerical  integration.  It  should  be 
noted  that  Eqs.  (1)  through  (3)  have  been  written  in  terms  of  the  dimensionless  size 
parameter  x  =  7rD/X,  where  X  is  the  wavelength  of  the  incident  radiation.  The  size 
parameter  is,  obviously,  the  ratio  of  the  particle  circumference  to  the  light  wavelength. 

The  differential  scattering  cross  sections  are  related  to  angular  Mie  intensity 
functions,  which  are  obtained  from  the  Mie  scattering  computer  codes,  according  to 


*Only  the  first  two  Stokes  parameters  of  scattered  light  —  the  intensity  and  the  degree-of- 
polarization  —  are  considered.  For  nonspherical  scatterers,  it  may  be  necessary  to  measure  the  remaining 
two  Stokes  parameters,  the  ellipticity  and  orientation  of  the  polarization  ellipse.  Such  additional 
measurements  are  also  desirable  when  dealing  with  multiple  scattering.  However,  in  this  work,  these 
additional  complications  are  ignored. 
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dtJ-j 

an 


(x,  9,  rj)  =  ij  (x,  e,  17)  Sin  2  ^/k2 


(4) 


and 


— f  (x,  -  ig  (*.  9,  17)  cos 2  ^/k2  (5) 

dU 

where  k  =  2tt/A  is  the  wave  number  of  the  incident  light.  In  these  equations,  7?  =  7jt  -  iic 
is  the  complex  index  of  refraction  of  the  particle  species  relative  to  that  of  the  medium 
in  which  they  are  embedded  (in  this  case,  air,  whose  index  is  unity). 

When  a  polydispersion  of  sizes  must  be  considered,  Eqs.  (1)  through  (3)  should  be 
replaced  by  analogous  forms  in  which  3n(x,  £)/9x  is  substituted  for  n(x,  fi),  and  the 
right-hand  sides  of  Eqs.  (1)  and  (2),  as  well  as  the  exponent  of  Eq.  (3),  are  then 
integrated  over  the  interval  o  <  x  <«.  Usually  the  scattering  volume  is  sufficiently  small 
that  variations  in  number  density  and/or  size  distribution  can  be  ignored  through  the 
distance  Ej  (the  length  of  the  focal  volume  corrected  for  viewing  angle  effects),  in  which 
case  3n(x,  £)/9x  can  be  replaced  in  Eqs.  (1)  and  (2)  by  nf(x),  where  f(x)  is  the 
normalized  size  distribution  function.  The  scattering  length,  Ej,  is  different  for  each  angle 
[i.e.,  Ej  =  E(flj)].  The  number  density,  n,  is  taken  to  be  the  average  of  the  particle  number 
density  throughout  the  scattering  volume. 

The  essence  of  the  AEDC  Mie  scattering  technique  is  to  measure  and  interpret 
suitably  defined  scattering  ratios  from  which  the  computer-based  analytical  procedures 
can  then  deconvolve  the  particle  size  distribution  and  determine  the  average  particle  size 

and  number  density.  The  use  of  such  ratios  avoids  the  difficulty  of  estimating  beforehand 

the  (usually  unknown)  particle  number  density  and  avoids  the  attendant  experimental 
difficulties  associated  with  making  accurate,  absolute  measurements  of  scattered  radiation 
signals. 

2.2  MEASUREMENT  OF  THE  AVERAGE  SIZE  PARAMETER 

USING  THE  MIE  SCATTERING  FUNCTIONS 

Several  functional  ratios  can  be  formed  that  yield  a  measure,  albeit  ambiguous,  of 
the  average  size  parameter  of  a  polydisperse  PSDF.  The  degree-of-polarization  function, 
an  example  of  such  a  ratio,  can  be  defined  as 

<  ij  (x,  0.,  Jjl  X  sin2  >- <ij  (x.  0f,  77)  >  <eoB2  > 

P  tf,.  *!.  X.  9)  -  <it  (?,(?.,,)>  <sin2  >+  <i2  (7,  <eos2  <|r,>  (6) 
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where 


and 


<ili2(*.  ffj,  Jj)> 


{*,  Ot,  7)  sin  8 


(7) 


(8) 


The  integration  limits  for  6  and  1//  are  defined  as 
Oi"  -  fli  +  tAfl/a) 
and 

^i"  =  v>j  +  (A(/r./2)  -  (A^/2) 

where  sin  6>j  A0t  A^j  represents  the  included  solid  angle  subtended  by  the  collection 
system  at  the  scattering  volume.  The  resulting  integrations  for  <sin2  ^j>,  <cos2  J/;>  are 
given  by 


^sin2  i/r^  =  sin2  ^r.  -  (cos  2i/ri/2&tfri)  (sin  Ai/r.  -  A^r.) 


and 


<cos2  v>j>  =  cos2  \b-t  +  (cos  2^./2A^j)  (sin  A -  A<£.)  (ga) 

Shown  in  Fig.  6  is  the  variation  of  P(0L,  x,  7)  with  size  parameter  x  for  water 
(H20)  at  two  backscattering  angles  near  the  rainbow  angle  for  H20.  The  azimuthal  angle 
of  7t/4  was  assumed. 

Comparison  of  Eqs.  (1),  (2),  (4),  and  (5)  with  Eq.  (6)  shows  that  the  experimentally 
determined  degree-of-polarization,  Pcxp  (01;  ^i)>  can  be  written  as 


P”P  {$.,  <(,.)  =  [<t>l  (*,,  v>j)  - 1»2  W,,  i  Wj.  0t)  +  »[•  *i>]  <9) 


If  a  value  of  PexP  is  assumed  (0j,  i/q),  for  which  the  parenthetical  parameters  x  and  rj 
have  been  suppressed,  it  can  be  seen  (Fig.  6)  that  ambiguity  results  in  the  determination 
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of  x.  For  example,  if  PexP  (0i5  i//j)  were  found  to  be  0.60  in  Fig,  6a,  the  value  of  the 
measured  size  parameter  would  be  5.15,  5.48,  6.2,  6.75,  7.48,  or  7.85.  Another  factor 
contributing  to  the  uncertainty  Ax  in  the  size  parameter  results  from  experimental 
uncertainties  in  the  measured  flux  values,  <3>],2-  Specifically,  the  uncertainty  5P  in 
pxp(0i,  J/i),  resulting  from  uncertainties  in  the  flux  values,  is  given  by 

<5P  (flj.  ^>/p  Wj,  ^  =  2  (B/A>[(5A/A)2  +  {SB/B)2]Vi/  [l  -  (B/A)2]  (10) 


where 


A  = 
and 

B  =  4*2 

Because  the  steepness  of  the  oscillatory  functions  shown  in  Fig.  6  is  greatly  dependent  on 
the  three  parameters  0,  x,  and  rj,  it  is  difficult  to  determine  the  exact  relationship 
between  the  imprecision  of  the  degree-of-polarization  measurements  and  the  propagated 
imprecision  of  the  inferred  average  size  parameter.  Nevertheless,  from  the  application 
described  in  Ref.  3,  it  was  found  that  an  imprecision  in  P  of  from  approximately  10  to 
15  percent  resulted  in  a  derived  size  parameter  imprecision  of  from  about  15  to  20 
percent. 

Because  of  the  non-uniqueness  demonstrated  in  Fig.  6  and  because  of  the  effect  of 
propagated  measurement  imprecision,  it  is  apparent  that  either  additional  angular 
measurements  of  the  degree  of  polarization  or  additional,  independent,  measured  ratios  of 
the  Mie  functions  are  required  to  decrease  the  ambiguity  in  x.  Two  such  functions 
(ANGER)  are  the  ratios  formed  by  the  polarization-specific,  angular  scattering  signal  and 
the  transmission  signal: 


1*1=  <iL  (x,  dit  i/frj,  q)>/x2  Qext  (x,  7))  (11) 

and 

R2  =  6i’  ^i>  7>>A2  (*>  7)  (12) 

where 

<*•  ?>A(D/2)2  =  Qscat  (x,  V)  +  Qabe  (3c,  ,)  (13) 
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is  the  extinction  efficiency.  The  corresponding  experimental  quantities  Re,xP  and  R§xp  are 
defined  to  be 


RfP  =  (ff/<8in2  ^)/< T)  (14) 

RrP  =  (^/<cos2^i>^i)(L/£i)[®2iei.^i)/<I.o]/(-fn  T)  (15) 


where  T  is  the  experimental  value  of  the  transmission  ratio  as  measured  for  an  effective 
extinction  length  L.  L  is  a  “lumped”  parameter  that  incorporates  the  effect  of  particle 
density  gradients;  it  is  defined  as 


L  = 


dn(x,  l) 
(?x 


— ?)  dxdf 

<?X 


Figures  7a  and  b  show  the  variations  of  Rt  and  R2,  respectively,  with  x  for  iron  in  a 
particular  scattering  geometry.  Again,  ambiguity  in  the  determination  of  x  results. 

Now,  however,  the  task  at  hand  is  to  determine  whether  common  x  values  result 
from  measurement  of  the  P®xp  and  R®*!  signals.  Prior  to  this  search  of  the  set  of 
allowed  x  values,  the  experimental  uncertainty  values  8R®X?  are  required  to  determine 
the  corresponding  uncertainty  in  x.  These  relations  can  be  shown  to  be 


~  -  [«»l,2/®L,2)2^St"T/'"T>2]'/' 

1 ,  £  L 


(16) 


where  §4>]  2  and  6SnT  are  the  estimated  imprecision  values  for  the  corresponding 
measured  quantities,  and  5L  is  the  uncertainty  in  L. 

The  combination  of  Pcxp  and  R“X2P  for  determining  x  was  used  in  Ref.  3;  the 
diameter  of  the  particles  measured  in  Ref.  3  was  found  to  be  0.8  ±0.1  p.  In  this  case, 
the  predominant  contribution  to  the  uncertainty  of  x  was  the  uncertainty  in  the 
extinction  length  L,  and,  relative  to  8L,  both  and  5T  were  negligible. 

For  those  situations  in  which  the  measurement  of  the  transmission  ratio  T  is 
impractical,  an  additional  scattering  function  ratio  can  be  defined  wherein  the 
transmission  measurement  is  replaced  by  a  backscattering  channel.  This 
Angular/Backscattering  ratio  (ANGBAR)  can  be  defined,  for  the  independent  polarization 
states,  as 
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Rfil  2  =  (x*  >?)  =  <>1,2  (x*  ®i'  17^/*2  (x«  ^  0?) 

where  Qg  is  the  backscattering  efficiency  factor.  The  experimental  counterparts  to  Eq. 
(17)  are  obvious. 

Once  the  average  size  has  been  determined  with  any  of  the  sizing  procedures 
discussed  here,  particle  number  density  can  be  determined  from  Eq.  (3),  if  transmission 
measurements  are  available.  The  result  is  a  spatially  averaged  number  density  along  the 
Line  of  sight.  Of  course,  an  alternative  is  to  measure  the  absolute  scattering  intensity.  If 
transmission  measurements  are  available,  then  the  spatially  averaged  number  density  of 
the  scatterers  can  be  determined  from  Eq.  (3)  by  convolution  with  the  cross  sections, 
such  as  those  shown  in  Fig.  8.  This  type  of  analysis  was  demonstrated  in  Ref.  3. 

It  can  be  seen  by  referring  to  Eqs.  (11),  (12),  and  (17)  that  the  ANGER  and 
ANGfiAR  scattering  ratios  depend  upon  the  extinction  and  backscattering  efficiency 
factors  Qexi  and  Qb,  respectively.  A  better  appreciation  of  these  scattering  ratios  is 
achieved  by  noting  some  of  the  qualitative  dependencies  of  Qext  and  Qb  on  the 
parameters  x  and  37.  Shown  in  Fig.  9  is  the  variation  of  Qext  with  x  for  transparent  (n  = 
0)  and  absorptive  (k  =£  0)  particles.  For  absorptive  particles  (Ref.  9)  Lim  Qext  =  2. 

^  *4  on 

Figure  9  shows  that  convergence  to  the  limit  is  rapid  for  x.  Clearly,  this  asymptotic 
convergence  is  the  mathematical  expression  of  the  transition  in  scattering  theory  from  the 
exact  Mie  theory  to  the  geometrical  optics  limit.  No  such  limit  occurs  for  transparent 
particles  (Ref.  12).  In  fact,  oscillations  about  the  limit,  2,  persist  for  size  parameters  as 
large  as  several  thousand.  Figure  9  shows  that  there  are  several  scales  of  oscillations 
corresponding  to  resonances  in  different  terms  in  the  Mie  series.  Surface  waves,  generated 
by  the  last  few  terms  in  the  Mie  series  just  before  convergence,  are  the  sources  of  the 
smallest  scale  of  oscillations.  These  oscillations  are  characterized  by  size  parameter 
periodicity  much  less  than  unity,  and  they  disappear  altogether  when  the  absorptive 
index  is  as  large  as  0.001  (Ref.  13). 

The  backscattering  efficiency  has  a  similar  geometrical  optics  limit: 

Lim  Qg  (x,  Tf)  =  |  (77 +  1)  |2  (18) 

X-*oc 

As  is  shown  by  Fig.  10,  this  limit  also  occurs  at  comparatively  small  size  parameters  for 
absorptive  particles  but  not  for  transparent  particles.  As  is  the  case  for  Qcxt,  the 
oscillations  become  damped  with  an  average  value  that  approximates  the  geometrical 
optics  limit  at  size  parameters  on  the  order  of  a  few  thousand  (Ref.  10).  For  a 
polydisperse  PSDF,  the  oscillations  persist  with  respect  to  the  modal  diameter,  as  is 
shown  by  Fig.  1 1. 
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As  a  consequence  of  this  behavior  of  the  efficiencies,  the  ANGER  and  ANGBAR 
functions  vary  considerably  for  transparent  and  absorptive  particles.  Figure  1  2  shows  the 
ANGER  functions  for  copper  particles  and  for  imaginary  particles  having  the  same  real 
refractive  index,  but  zero  absorptive  index.  These  results  indicate  that,  because  of  the 
sensitivity  of  this  function  to  the  absorptive  index,  Mie  scattering  can  be  used  to 
determine  7j  in  an  iterative  fashion. 

Unlike  Qex  t  and  QBi  the  degree-of-polarization  measurements  are  insensitive  to  the 
absorptive  index,  but  they  are  quite  sensitive  to  the  real  refractive  index.  This  can  be 
seen  by  comparing  Fig.  6  with  Figs.  13  and  14,  which  show  the  degree-of-polarization 
versus  size  parameter  for  copper  and  iron,  respectively,  at  the  same  conditions  as  for  Fig. 
6.  Again,  these  results  suggest  that  the  real  refractive  index  can  be  determined,  iteratively, 
from  analysis  of  Mie  scattering  measurements.  For  iterative  determination  of  the 
refractive  indices,  angular  scattering  measurements  should  be  concentrated  near  the 
rainbow  angle  of  the  suspected  particle  constituents.  The  rainbow  angle  is  a  function  of 
only  the  real  refractive  index  in  the  geometrical  optics  limit.  As  the  size  parameter 
decreases,  the  rainbow  angle  increases  slowly  until  it  ceases  to  have  physical  significance. 
This  occurs  when  the  particle  diameter  becomes  too  small  to  permit  internal  reflections 
to  occur  (i.e.,  the  radiation  wavelength  either  equals  or  exceeds  the  particle  radius). 

The  preceding  discussion  of  properties  of  the  Mie  scattering  ratios  used  in  the  AEDC 
particle  diagnostics  scheme  was  illustrated  by  the  use  of  calculations  that  assume  an 
azimuthal  angle  of  jt/4  and  solid  angles  equal  to  those  stated  in  Ref.  3.  However,  when  it 
is  suspected  that  only  submicron  particles  are  present,  then  degree-of-polarization 
measurements  made  at  the  scattering  and  azimuthal  angles  of  irjl  and  either  0  or  ir, 
respectively,  can  be  extremely  useful,  as  shown  by  Fig.  15.  These  graphs  illustrate  the 
transition  from  Rayleigh  to  Mie  scattering;  this  behavior  permits  unique  determination  of 
the  average  size  of  particles  for  size  parameters  between  limits  that  depend  on  the 
collection  optics  solid  angles  and  the  real  refractive  index.  Figure  1 5  shows  the  effect  of 
progressively  increasing  the  subtended  solid  angle.  The  basis  for  this  method  of  size 
determination  can  be  explained  as  follows:  In  the  Rayleigh  scattering  limit,  the  parallel 
polarization  state  is  completely  suppressed  at  a  scattering  angle  of  ir/2,  as  shown  in  Figs. 
16  and  17  (Ref.  7).  As  shown  in  Fig.  16,  the  parallel  Mie  intensity  decreases  very  rapidly 
as  the  size  parameter  approaches  zero.  In  contrast,  as  shown  in  Fig.  1 7,  the  perpendicular 
Mie  function  decreases  much  less  rapidly  with  size  parameter.  The  use  of  incident 
polarization  parallel  or  anti-parallel  to  the  scattering  plane,  however,  suppresses  the 
perpendicular  component,  as  suggested  by  Eq.  (4),  except  for  terms  proportional  to  the 
subtended  azimuthal  angle.  [See  Eq.  (8a).]  With  a  decrease  in  particle  size,  the  degree  of 
polarization  increases  from  its  asymptotic  limit  of  -1  because  the  parallel  polarization 
component  becomes  comparable  to  the  azimuthally  suppressed  perpendicular  component. 


16 


AEDC-TR-80-3 


I 

This  behavior,  shown  in  Fig.  1 8,  results  from  the  transition  of  the  scattering  process  into 
the  Rayleigh  scattering  regime.  Reference  7  describes  the  use  of  this  transition  to 
determine  the  size  of  very  small  droplets  present  in  the  exhaust  of  a  liquid-propellant 
rocket  engine.  In  this  configuration,  measured  degree-of-polarization  values  greater  than  -1 
are  indisputable  evidence  that  the  PSDF  occupies  predominantly  the  size  range  below  the 
limit  at  which  asymptotic  behavior  begins. 

3.0  SIZE  DECONVOLUTION  ANALYSIS 

A  histogram  representation  of  the  PSDF  can  be  obtained  with  use  of  multiangular 
measurements  of  the  Mie-scattered  signals.  The  PSDF  deconvolution  of  the  integral 
scattering  equations  Jthe  polydisperse  forms  of  Eqs.  (1)  and  (2)]  can  be  accomplished,  in 
principle,  with  a  direct  inversion  of  the  equations  as  well  as  with  a  multivariate  nonlinear 
regression  analysis.  These  options  are  shown  in  Fig.  4. 

The  inversion  scheme  discussed  in  Appendix  A  yields  a  histogram  representation 
of  the  PSDF,  f(x).  This  representation  is  defined  for  each  of  several  equal-sized  bins  as 
the  result  of  integrating  the  PSDF  over  the  size  bin.  Thus,  the  results  of  deconvolution 
represent  the  mole  fraction  of  particles  in  each  size  bin.  Consequently,  the  histogram 
approximation  to  the  PSDF  is  defined  as  the  result  of  dividing  the  deconvolved  mole 
fractions  for  each  size  bin  by  the  bin  size.  The  DECRATS  and  the  deconvolved  histogram 
are  shown  and  discussed  in  Appendix  A  for  this  deconvolution  using  simulated  scattering 
data.  The  ultimate  size  resolution  (the  number  of  size  bins  in  the  histogram)  that  is 
attainable  by  this  size  deconvolution  procedure  is  limited  by  the  level  of  both 
experimental  and  computational  error  in  each  DECRAT;  this  is  detailed  in  Appendix  A. 

One  might  generally  expect  accurate  deconvolution  into  a  number  of  size  bins  equal 
to  the  number  of  forward  scattering  angles  plus  twice  the  number  of  backscattering 
angles  (excluding  the  reference  angle)  at  which  scattering  data  are  taken.  Actually,  the 
number  of  bins  that  can  be  deconvolved  accurately  is  somewhat  smaller  than  this 
approximation.  Whether  an  accurate  deconvolution  into  a  specified  number  of  size  bins 
can  be  accomplished,  for  a  given  experimental  configuration,  can  be  determined,  as 
explained  in  Appendix  A,  by  comparing  eigenvalues  of  the  Mie-scattering  kernels  (the 
results  of  quadrature  of  the  Mie  functions  in  each  size  bin)  with  the  computed  residual 
ereor  which  is  obtained  without  smoothing  the  kernels.  If  this  test  fails,  then  the  number 
of  size  bins  must  be  reduced. 

The  choice  to  invert  the  scattering  equations  by  this  deconvolution  procedure,  rather 
than  to  invert  directly  by  nonlinear  regression,  was  made  on  the  basis  of  the  strong 
non-monotonicity  of  the  Mie  functions.  For  such  functions  unique  results  are  difficult  to 
obtain  by  nonlinear  regression  analysis. 
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Another  technique  is  being  developed  that  will  combine  the  virtues  of  both  the 
Fredholm  deconvolution  and  direct  regression  procedures.  An  existing  AEDC  multivariate 
nonlinear  regression  scheme  (Ref.  14)  is  being  adapted  for  these  purposes.  The  resulting 
overall  regression  can  be  organized  into  the  following  steps: 

1.  The  histogram  approximation  to  the  PSDF  is  obtained,  as  described  in 
Appendix  A. 

2.  Parameters  of  the  PSDF  are  obtained  by  regression  of  the  histogram 
approximation. 

3.  The  Mie  scattering  functions  corresponding  to  the  regressed  PSDF  are  then 
computed  and  compared  with  the  experimental  data  to  see  whether  the 
computed  scattering  ratios  lie  within  the  experimental  imprecision  bands. 

4.  Numerical  derivatives  can  also  be  evaluated  for  the  direct  regression  of  the 
scattering  equations,  and  direct  regression  of  the  scattering  data  can  be 
undertaken,  using  the  results  of  the  previous  regression  as  initial  guesses  to 
the  full-scale,  direct  regression. 

Obviously,  steps  3  and  4  must  be  repeated  until  a  satisfactory  regression  is  obtained. 
However,  the  stability  characteristics  of  this  technique  are  at  present  unknown  and  must 
be  determined  in  the  course  of  developing  this  method  of  inversion. 

4.0  DETERMINING  THE  AVERAGE  SIZE  PARAMETERS  AND 
MOLE  FRACTIONS  OF  PARTICULATE  MIXTURES 

In  principle,  it  is  possible  to  use  Mie  scattering  for  determining  the  average  particle 
size  parameters  and  species  mole  fractions  of  a  polydisperse,  multispecies  particulate 
mixture  of  known  indices  of  refraction.  The  only  requirement  is  that  there  be  an 
adequate  number  of  Mie  scattering  measurements  consisting  of  the  polarization  state, 
angular  distribution,  and  wavelength  variation  of  the  scattered  radiation.  These 
independent  variables  may  be  taken  singly  (such  as  taking  only  the  angular  variation 
measurements)  or  in  combinations  (such  as  taking  both  wavelength  and  polarization 
variations  of  the  angular  scattering  dependence  of  the  Mie  signal).  For  simplicity,  the 
following  derivation  will  assume  the  use  of  both  states  of  polarization,  multiangular 
measurements,  and  the  use  of  a  single  wavelength.  In  consideration  of  the  existing 
uncertainties  in  deconvolving  the  PSDF  of  a  single-species  particulate  sample,  the  initial 
development  of  a  multispecics  deconvolution  has  been  the  application  of  Mie  scattering 
theory  to  a  binary  particulate  mixture  of  known  refractive  indices;  the  individual  species 
distribution  functions  of  this  mixture  are  treated  as  though  they  were  monodisperse  with 
respect  to  size  parameter. 
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One  can  easily  generalize  the  degree  of  polarization,  P,  as  given  by  Eq.  (6),  to  apply 
to  a  binary  particulate  mixture  of  mole  fractions  and  size  parameters,  as  and  xs,  s  =  1.2, 
respectively.  For  scattering  angles  0j,  ~  *14,  the  degree  of  polarization  can  be  written 

as 

P(0j,it/4)=  2[<I1S-Ii>y£[ais  +  J,K]  (19a) 

8=1 

S-KS  'V 

where  the  elements  of  the  matrices  I  and  J  are  given  by 

=  <‘1(6ir>?s,  x5)>  (19b) 

and 

Jis  =  ^'2  ’Jg' xs)  ^  (19c) 

Obviously,  the  range  of  the  subscript  i  depends  upon  the  number  of  scattering  angles 
used  for  the  measurement.  If  P(04,  tt/4)  is  defined  as  Pj,  then  Eq.  (19a)  can  be  rewritten 
as 


2 


S=1 


0 


(20a) 


where 


t„  -  ffi-inl.  +  ip.  +  in,.  wow 

When  two  angular  measurements  have  been  obtained  (i  -  1,2  for  the  binary 
mixture),  then  the  set  of  homogeneous  equations  of  Eq.  (20a),  usually  an  overdetermined 
set,  becomes  a  completely  determined  set.  By  expanding  Eq.  (20a)  for  the  stated  case, 
one  finds  that 


Tllal  +  T12  °2 


0 


(20c) 


and 


T21  al  +  T22  °2  “  0 


(20d) 
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The  ratio  of  mole  fractions  0  =  ai  ja-i  is,  of  course,  found  to  be 

fi  ■  -Tij/t,,  -  -T22/r„  (21a) 

It  is  appropriate  to  recall  that  Pj  is  an  experimentally  determined  parameter.  With  use  of 
Eq.  (20b),  one  sees  that  (3  depends  not  only  on  Pj  but  also  on  the  quantities  IjS  and  JjS, 
both  of  which  are  functions  of  the  unknown  size  parameters  xs,  s  =  1,2.  Consequently, 
an  iterative  solution  of  Eq.  (21a)  for  0  is  required  for  which  values  of  the  size  parameter 
pair  Xu  are  assumed.  Eqs.  (19b)  and  (19c)  are  used  to  compute  IjS  and  JjS,  respectively. 
Finally,  since  the  solutions  for  as,  s  -  1,2,  are  mole  fractions  of  the  binary  mixture,  the 
solution  set  {a,}  is  subject  to  the  following  constraints: 

°s  >  o  (21b) 


and 


S  «  =1  (21c) 

S  =  1  S 

Since  the  process  of  solving  Eq.  (20a)  is  both  numerical  and  iterative,  the  parameters 
jx^k ^ ,  s  =  1,2,  are  defined  to  be  the  kth  set  of  assumed  size  parameters,  and  xs,  s  =  1,2, 
represent  the  actual  size  parameters  of  the  mixture.  Further,  in  anticipation  of  the 
round-off  error  that  will  be  encountered  in  the  numerical  solution,  the  constraint  of  Eq. 
(21c)  is  relaxed  to  be 


2 
2 
s -1 


<<  1 


where  the  superscript  k  on  all  quantities  denotes  their  values  appropriate  to  the  kth 
choice  of  the  set  {xfJO}. 

The  iterative  determination  of  the  actual  size  parameters  of  the  mixture  is  begun  by 
casting  Eq.  (21a)  into  the  form  of  an  inequality 


P  _p<k>  M(k>  P  -P(k> 
l  11  1  l  12 


P 


t 


P<t) 

22 


V|(k)  p  _  p(k) 
2  2  22 


<£  <<1 


(22) 
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where  the  auxiliary  functions  >  and  P0j>  are  defined  for  the  kth  set  of  size  parameters, 
{xW},  as 


If  the  kth  set  of  size  parameters,  satisfies  Eq.  (22),  the  corresponding  species 

mole  fractions  are  found  from  Eq.  (21a)  and 


and  =  - - - 

l  +  0'l> 


(24) 


a  early,  the  correct  size  parameters  and  species  mole  fractions  are  those  choices  that 
satisfy  both  Eqs.  (21  b)  and  (22). 

Because  of  uncertainties  in  the  measurements,  there  will  be  a  range  of  size 
parameters  and  species  mole  fractions  that  satisfies  both  Eqs.  (21b)  and  (22).  The 
propagation  of  data  imprecision  into  the  analysis  is  taken  into  account  when  choosing  the 
set  of  experimental  values  { Pj }  so  that  they  lie  within  the  imprecision  intervals  stated  in 
Eq.  (10).  If  the  jth  such  set  is  denoted  then  each  P(p  mus^  satisfy 


PexP  -  BP,  <  <  P?xp  +  5P< 

1  —  1  —  1  i 


For  each chosen  the  entire  algorithm  must  be  used  to  determine  the  allowed  particle 
sizes  and  species  mole  fractions.  The  number  of  samples  of  {P^}  which  should  be 
investigated  represents  a  compromise  between  (1)  the  requirement  of  a  large  number  of 
samples  to  assure  adequate  sampling  statistics  and  (2)  the  limitations  of  the  computer. 

This  analysis  was  applied  in  Ref.  3  (see  Appendix  B)  to  determine  species 
concentrations  of  copper  particles  and  tungsten  particles  in  the  AEDC  von  Karman  Gas 
Dynamics  Facility  (VKF)  Tunnel  F.  The  analysis  indicated  that  the  Mie  scattering 
measurements  performed  at  two  angular  stations  were  consistent  with  particulate  mole 
fraction  ranges  of  from  100%  copper/0%  tungsten  to  79%  copper/21%  tungsten.  Posttest 
analysis  of  a  sample  of  collected  particles  indicated  the  composition  to  be  75% 
copper/25%  tungsten.  In  addition,  the  two-species  average  size  parameters  were  within  the 
uncertainty  range  of  the  average  size  parameter  for  the  case  of  100%  copper,  and  the 
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total  mass  of  particulate  matter,  estimated  with  use  of  the  analytical  procedures  reported 
herein  was  consistent  with  measurements  of  electrode  metal  loss  in  the  Tunnel  F 
arc-chamber,  the  source  of  the  particles. 

5.0  CONCLUDING  REMARKS 

This  report  has  described  and  summarized  the  systematic  development  of  the 
analytical  procedures  used  for  Mie  scattering  particle  diagnostics.  It  is  believed  that  this 
program  will  prove  to  be  useful  in  a  variety  of  applications  at  AEDC  and  other  agencies. 
The  strengths  of  the  approach  are  (1)  its  ability  to  resolve  the  PSDF  for  an  ensemble  of 
many  particles  present  in  the  focal  volume  of  each  detector,  (2)  time-resolution  capability 
on  the  order  of  miliseconds,  (3)  ability  to  resolve  particle  chemical  species  distribution  (if 
species  refractive  indices  are  known),  (4)  limited  capability  to  determine  effective  indices 
of  refraction  in  an  iterative  fashion  for  samples  of  unknown  properties,  and  (5)  the 
ability  for  implementation  using  available  laser  sources,  detectors,  and  associated 
electronic  and  data  recording  systems.  The  program  relies  heavily  on  analytical  procedures 
performed  with  a  large-scale  computer,  but  aspects  of  it  have  also  been  successfully 
performed  on  a  microcomputer,  the  results  of  which  demonstrate  that  the  overall 
diagnostics  scheme  is  versatile  enough  to  permit  incorporating  portions  of  it  into  portable 
instrumentation  and  data  analysis  packages. 

Future  work  will  concentrate  on  improvement  and  additional  testing  of  the 
techniques  developed  and  completion  of  those  portions  of  the  program  that  remain  under 
development.  Incorporation  of  multispectral  capability  into  the  basic  diagnostic 
procedures  discussed  here  will  also  be  investigated.  Other  potential  extensions  of  the 
program  include  such  obvious  areas  as  Raman  scattering  and  fluorescence  from  structured 
particles.  In  addition,  generalization  to  permit  diagnostics  of  aspherical  particles  is 
currently  under  consideration,  and  the  effectiveness  of  these  basic  diagnostic  methods  in 
multiple-scattering  environments  will  be  considered. 
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Figure  1 .  Comparison  of  Mie  theory  with  experimental  scattered 

light  intensity  from  monodisperse  polystyrene  latex  spheres 
in  water  suspension  (from  Ref.  4). 
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Artie  Scattering  Intensity  Function  for  Perpendicular 
Polarization  State,  dimensionless 


Scattering  Angle,  deg  Scattering  Angle,  deg 

a.  Perpendicular  polarization  state  b.  Parallel  polarization  state 


Figure  2.  Comparison  of  Mie  theory  with  experimental  scattered  light  intensity 
from  a  narrow  polydispersion  of  spherical  aerosols. 
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Figure  5.  Spherical  coordinate  geometry  for  the  multiangular 
scattering  of  light  by  a  dielectric  sphere. 
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Size  Parameter,  dimensionless 

b.  6  =  161.90  deg 

Figure  6.  Size  dependence  of  the  degree  of  polarization  of  light 
scattered  by  water  droplets  in  the  geometry  of  Ref.  3. 
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ANGER  Functions  tor  Parallel  Polarization  ANGER  Functions  for  Perpendicular  Polarization 

State,  dimensionless  State,  dimensionless 
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Size  Parameter,  dimensionless 

a.  Perpendicular  polarization 


Size  Parameter,  dimensionless 


b.  Parallel  polarization 

Figure  7.  Size  dependence  of  ANGER  functions  for  iron  In  the 
geometry  of  Ref.  3. 
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a.  Carbon  soot  —  v  ~  1-56  —  0.46  i 


Size  Parameter,  dimensionless 


b.  Calclte  (CaCOs)  —  t\  -  1 .6  —  0  i 

Figure  9.  Size  dependence  of  the  extinction  efficiencies  for 
carbon  soot  and  calcite:  visible  light. 
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Size  Parameter,  dimensionless 

a.  Carbon  soot  -  rj  =  1 .56  -  0.46  i 


20 . 0  «— 


16,0  h 


12.0  I- 


0.3 


12.3  24.3  36.3 

Size  Parameter,  dimensionless 


48.3 


b.  Calcita  (CaC03)  -  rj  =  1 .6  -  0  i 
Figure  10.  Size  dependence  of  the  backscattering  efficiencies 
for  carbon  soot  and  calcite:  visible  light. 


Note:  isometric  Standard  Deviation  -  2.0, 
Log-Normal  Size  Distribution 


I  t  I  I 

1.0  3.Q  3.0  4.0 

Kodal  Diameter,  u. 


Figure  1 1 .  Dependence  of  backscattering  efficiency  on  modal  diameter 

for  polydisperse  soot  particles  at  visible  and  infrared  wavelengths 


ANGER  Functions  far  Perpendicular  Polarization  ANGER  Functions  tor  Perpendicular  Polarization 

State,  dimensionless  State,  dimensionless 
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Size  Parameter,  dimensionless 


b.  T t  =  1.06  -  2.31  i 

Figure  1 2.  Effect  of  absorptive  index  on  size  dependence  of  ANGER  functions 
for  copper  (perpendicular  polarization)  in  the  geometry  of  Ref.  3. 
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Size  Parameter,  dimensionless 

a.  8  =  146.50  deg 


Size  Parameter,  dimensionless 


b.  6  =  161 .90  deg 

Figure  1 3.  Size  dependence  of  the  degree  of  polarization  of  light  scattered 
by  copper  spheres  in  the  geometry  of  Ref.  3. 
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Site  Parameter.  dimcnsion.pw  Size  Para-.eter.  iimensnmless 

a.  A4r  =  1  deg,  AO  =  1  deg  b.  Aty  =  2  deg,  AS  =  2  deg 


Size  Parameter,  dimensionless  Size  Parameter,  dimensionless 

c.  AV  =  3  deg,  AO  =  3  deg  d.  A'l'  =  4  deg,  AO  =  4  deg 

Figure  1 5.  Effect  of  collection  optics  subtended  solid  angle  on  size  dependence  of 
the  degree  of  polarization  for  light  scattering  from  AI2O3  spheres  at 
90-deg  scattering  angle  and  incident  polarization  parallel  to  the 
scattering  plane. 
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Log  10  of  Hie  Scattering  Functions 
for  Parallel  Polarization  State 


Figure  16.  Size  dependence  of  the  parallel  Mie  scattering  Intensity 
functions  in  the  geometry  of  Ref.  7. 
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Log  10  of  Mie  Scattering  Functions 
'or  Perpendicular  Polarization  State 


Size  Parameter,  dimensionless 

Figure  17.  Size  dependence  of  the  perpendicular  Mie  scattering  intensity 
functions  in  the  geometry  of  Ref.  7. 
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Degree  of  Polarization,  dimensionless 


Size  Parameter,  dimensionless 

Figure  1 8.  Size  dependence  of  the  degree  of  polarization  of  light  scattered 
by  very  small  droplets  in  the  geometry  of  Ref.  7. 
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APPENDIX  A 


SIZFRED: 

A  COMPUTER  PROGRAM  FOR  SIZE  DISTRIBUTION 
DECONVOLUTION  BY 

FREDHOLM  INVERSION  OF  MIE  SCATTERING  DATA 


B.  P.  Curry 
ARO,  Inc. 
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SUMMARY 

This  appendix  is  a  user’s  manual  for  code  SIZFRED,  a  program  developed  at  AEDC 
for  deconvolution  of  particle  size  distribution  histograms  obtained  by  single-wavelength, 
multiple-angle  Mie  scattering  measurements.  The  experimental  inputs  to  the  code  are 
Deconvolution  Ratios  (DECRATS),  which  are  defined  to  be  the  calibrated  signal  in  each 
orthogonal  polarization  state  at  each  angular  position,  normalized  by  the  unpolarized 
scattering  signal  in  the  90-deg  scattering  channel.  The  incident  laser  polarization  is 
assumed  to  be  45  deg  with  respect  to  the  scattering  plane. 

Use  of  this  code  permits  particle  size  deconvolution  into  as  many  as  30  size  bins  by 
using  a  maximum  of  40  input  ratios  (DECRATS).  The  code  incorporates  error 
propagation  analysis,  and  it  includes  a  provision  for  estimating  the  upper  bounds  on  the 
imprecision  of  each  separate  size  bin  in  the  deconvolved  histograms.  This  option  is 
activated  whenever  the  estimated  imprecision  of  each  DECRAT  is  punched  on  the  same 
input  card  as  the  corresponding  DECRAT.  The  resulting  deconvolution  is  a  least-squares 
solution  in  which  greatest  weight  is  given  the  DECRATS  having  least  imprecision.  If 
DECRAT  imprecision  values  are  not  specified,  the  code  defaults  to  an  option  in  which  all 
DECRATS  are  qually  weighted. 


A1.0  INTRODUCTION 

Development  and  application  of  rapid,  time-resolved  in-situ  particle  diagnostics  are 
essential  for  a  complete  understanding  of  rocket  exhausts,  wind  tunnel  flows,  and  other 
particulate-laden  flows.  Accurate  diagnostics  of  polydisperse  particle  size  distributions 
with  diameters  mainly  in  the  micron  and  submicron  range  requires  the  use  of 
experimental  Mie  scattering  measurement  techniques  and  appropriate  computational 
procedures  to  deconvolve  the  distribution.  This  study  was  undertaken  to  provide  a  size 
deconvolution  code  (SIZFRED)  that  incorporates  error  propagation  analysis  and  includes 
a  procedure  for  identifying  the  optimum  smoothing  level  of  the  deconvolved  particle  size 
histogram. 

This  appendix  explains  the  theory  and  application  of  code  SIZFRED  to 
deconvolving  particle  size  distribution  functions  from  angular  scattering  data  taken  in 
both  orthogonal  polarization  states.  These  states  are  defined  as  the  scattered-light  states 
in  which  the  polarization  vectors  are  perpendicular  to  and  parallel  to  the  plane  containing 
the  incident  and  scattered-light  propagation  direction  (i.e.,  the  scattering  plane). 
SIZFRED  assumes  that  the  data  in  each  polarization  state  from  each  angular  channel 
have  been  normalized  by  the  unpolarized  signal  from  the  90-deg  scattering  channel. 
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The  deconvolution  ratios  are  input  as  follows:  First,  ratios  for  the  perpendicular 
state  are  read  into  elements  (1,  N)  of  the  DECRAT  vector;  the  remaining  elements  (N  + 

1 ,2N)  are  the  ratios  corresponding  to  the  parallel  polarization  state.  Thus,  2N  DECRAT 
cards  must  be  read  in.  (N  is  the  number  of  angular  channels,  excluding  the  90-deg 
reference  channel.) 

This  program  differs  from  an  earlier  computation  in  three  important  respects: 

1.  It  has  a  provision  for  including  estimated  imprecision  in  the  input 
deconvolution  ratios. 

2.  It  weights  the  scattering  kernel  matrix  elements  with  the  reciprocal  of  the 
variance  of  the  appropriate  signal,  and  it  determines,  with  standard 
statistical  procedures,  the  propagation  of  these  estimated  imprecision 
values  into  the  final  solution  imprecision  bands.  Thus,  1 -sigma  input 
imprecision  values  produce  solution  imprecision  bands  which  are  the  upper 
bounds  to  the  corresponding  68-percent  confidence  limits. 

3.  It  provides  for  estimating  the  optimum  level  of  smoothing  through 
locating  the  minimum  in  either  or  both  of  two  goodness-of-fit  criteria 
(GOOF).  If  there  is  no  minimum  in  either  parameter,  then  the  GOOF 
parameters  decline  rapidly  at  first  and  then  reach  an  asymptotic  condition 
as  a  function  of  the  smoothing  parameter.  The  apparent  inflection  point  at 
which  this  asymptotic  behavior  occurs  can  be  considered  the  optimum 
smoothing  point,  provided  that  it  occurs  in  a  region  in  which  the  solution 
is  non-negative. 

In  earlier  simulation  calculations  with  similar  codes,  the  optimum  smoothing  point 
nearly  always  has  occurred  at  a  smoothing  level  such  that  all  parameters  of  interest 
(including  the  solution  vector)  vary  slowly  with  the  smoothing  parameter.  Thus,  the 
necessity  for  subjective  judgments  regarding  when  to  stop  smoothing  has  been  somewhat 
alleviated  in  this  code,  as  compared  to  its  predecessors. 

This  appendix  explains  the  use  of  the  code  and  gives  several  examples  of  successful 
deconvolutions,  some  of  which  include  the  insertion  of  simulated  “experimental  error." 

A2.0  GENERAL  THEORY 
A2.1  LIGHT-SCATTERING  RATIOS 

For  a  single-species,  polydisperse  ensemble  of  dielectric  spheres,  the  ratio  of  power 
scattered  from  a  localized  focal  volume  (determined  by  the  collection  optics  of  the 
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experiment)  to  the  incident  power  can  be  written  (Refs.  A-l  through  A-4)  in  each 
orthogonal  polarization  state  as  follows: 


0 


^sin2  >/.  f(x)<ij  (x,  0a,jj)>dx 


^2^8’  4>sr  y) 


f(x)<i2  (x,  0S, 


(A-l  a) 


(A-lb) 


where  n  is  the  number  density  of  the  particles,  £  is  the  projection  of  the  focal  volume  on 
the  incident  laser  axis  (scattering  length),  k  is  the  incident  light  wave  number  (inverse 
centimeters),  f(x)  is  the  normalized  size  distribution  function,  and  x  is  the  size  parameter 
(ratio  of  circumference  to  wavelength)  for  a  sphere  of  diameter  D.  The  angles  9  and  i^, 
defined  in  the  context  of  Fig.  A-l,  are,  respectively,  the  polar  scattering  angle  and  the 
azimuthal  angle  for  a  particle  that  is  located  at  the  origin  and  illuminated  from  below 
(along  the  z-axis)  by  a  laser  whose  electric  field  is  polarized  along  the  x-axis.  r?  is  the 
particle’s  complex  index  of  refraction. 

For  purposes  of  analysis,  computer-generated  Mie  scattering  functions,  integrated 
over  the  appropriate  subtended  angular  range,  are  defined  by  the  relations 
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Although  the  above  analysis  is,  strictly  speaking,  only  appropriate  for  spheres,  various 
modifications  to  Mie  scattering  theory  have  been  developed  to  simulate  the  scattering  by 
randomly  oriented,  polydisperse  nonspherical  scatterers.  In  such  cases,  the  appropriately 
modified  Mie  functions  correspond  to  the  use  of  surface-equivalent  diameters. 

Because  of  the  occurrence  of  the  (unknown)  number  density  in  Eq.  (A-l)  and  the 
difficulty  of  obtaining  accurate  absolute  intensity  measurements,  deconvolutions  are 
performed  by  using  ratios  of  experimental  scattering  signals.  Specifically,  the  data  inputs 
are  deconvolution  ratios  (DECRATS)  which  are  defined  as  the  ratio  of  the  scattering 
signal  in  a  given  angular  channel  in  each  separate  polarization  state  to  the  unpolarized 
90-deg  scattering  signal.  The  90-deg  signal  was  chosen  as  a  reference  denominator  because 
of  its  comparative  insensitivity  to  the  particle  distribution  function  (PSDF)  of  the  sample 
size.  The  implications  of  this  form  of  DECRAT  are  explored  in  Ref.  A-4. 

A2.2  CONSTRAINED  LINEAR  INVERSION 

Insofar  as  the  size  dependence  of  the  reference  denominator  can  be  ignored  (thus 
treating  the  90-deg  signal  as  a  normalization  constant),  Eq.  (A-l)  has  the  form,  with 
respect  to  the  size  parameters,  of  Fredholm  equations  of  the  first  kind.  Phillips  (Ref. 
A-5)  and  Twomey  (Refs.  A-6  through  A-8)  have  developed  a  basic  procedure  for  inverting 
such  equations,  taking  full  account  of  computational  errors  and  random  experimental 
“noise”  in  the  measured  inputs  (DECRATS).  Without  allowance  for  these  errors,  the 
PSDF  solutions  of  Eq.  (A-l)  are  unstable  and  exhibit  physically  unrealistic  oscillations 
and  negative  values.  If  the  vector  of  input  DECRAT  values  is  denoted  as  G,  then  Eq. 
(A-l)  can  be  rewritten,  allowing  for  error,  in  matrix  form  as  (using  the  Einstein 
summation  convention  for  repeated  subscripts) 

Qki(Gk  +  ck)  =  Gi  =  Cfri),  Fj  -  p  f(x)dx  (A-3) 

J  xj  -l 

In  Eq.  (A-3),  M  is  the  symmetric  scattering  kernel  matrix,  defined  as 

Mi,  -  0,j  -  f  Kfx.Xitdx  (A-4) 

"j-1 

The  vector  E,  introduced  in  Eq.(A-3),  is  comprised  of  the  residual  error  associated  with 
each  DECRAT.  These  errors  are  computed  from  Eq.  (A-3)  after  the  solution  vector  F  has 
been  computed.  The  norm  of  the  residuals  is  also  computed. 
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In  Eq.  (A-4),  variables  y  and  x  represent,  respectively,  the  angular  variables  and  the 
size  parameter.  The  specific  relations  that  relate  these  equations  to  the  light-scattering 
deconvolution  problem  are  the  following: 
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where  N  is  the  number  of  angular  channels.  The  associated  scattering  kernels  are  defined 
as 
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For  practical  reasons,  the  semi-infinite  limits  on  the  integrals  in  Eq.  (A-l)  have  been 
replaced  by  minimum  and  maximum  size  parameters  in  Eqs.  (A-7)  and  (A-8).  Although 
these  equations  include  general  azimuthal  forms,  it  is  assumed  henceforth  that  the 
incident  laser  polarization  axis  is  oriented  at  45  deg  with  respect  to  the  scattering  plane 
and  that  the  subtended  azimuths  (but  not  subtended  scattering  angles)  are  the  same  for 
all  angular  channels.  If  such  is  not  the  case,  the  deconvolution  code  must  be  modified 
accordingly. 

Twomey  (Ref.  A-5)  has  shown  that  the  constrained  least-squares  solution  of  Eq. 
(A-3)  can  be  stated  for  a  given  level  of  mean-square  residual  error  as 


Fi  =  K*1  W 

Rik  =  Mik  +  yHit 


(A-9) 


The  function  of  the  constraint  is  to  select  one  of  the  infinite  manifold  of  solutions  which 
are  possible  as  a  result  of  the  occurrence  of  error  in  Eq.  (A-3).  Here,  7  is  a  Lagrange 
multiplier  (the  smoothing  parameter);  the  matrix  H  is  determined  by  the  form  of  the 
constraint  chosen,  which,  in  this  case,  is  maximum  smoothness  of  the  solution  and  is 
internally  generated  in  the  code.  The  matrices  Q,  M,  and  R  represent,  respectively,  the 
unsymmetric  scattering  kernel,  symmetric  scattering  kernel,  and  smoothed  symmetric 
scattering  kernel  matrices. 

Twomey  originally  intended  that  the  optimum  value  of  the  smoothing  parameter  be 
determined  by  increasing  it  from  zero  until  the  computed  mean-square  error  equaled  its 
corresponding  experimental  value.  It  turns  out  that  the  computed  error  is  not  monotonic 
in  its  variation  with  respect  to  the  smoothing  parameter  (Ref.  A-9),  and  simulated 
deconvolutions  with  included  “noise”  in  the  inputs  show  that  the  optimum  smoothing 
parameter  corresponds  to  a  level  of  computed  mean-square  error  that  is  considerably 
smaller  than  the  simulated  “experimental”  noise  input.  Nevertheless,  when  a  known 
distribution  function  is  deconvolved  from  computed  scattering  data,  a  distinct  minimum 
in  the  variance  of  the  distribution  is  found  (Ref.  A-9).  It  is  fortunate  that  optimum 
smoothing  generally  occurs  in  a  region  in  which  the  variation  of  the  deconvolution 
solution  with  7  is  relatively  slow.  An  approximate  procedure  for  locating  the  optimum 
smoothing  level  will  be  discussed  later. 

The  need  for  smoothing  arises  because  of  the  sensitivity  of  Eq.  (A-3)  to  error. 
Twomey  (Ref.  A-7)  has  shown,  that  whenever  any  eigenvalue  of  the  symmetric  kernel 
matrix  is  not  significantly  larger  than  the  ratio  of  the  norm  of  the  computed  error  to  the 
norm  of  the  solution  vector,  then  the  kernel  is,  for  practical  purposes,  singular,  and  the 
number  of  independent  measurements  must  be  reduced  by  one  for  each  eigenvalue  for 
which  this  is  true.  The  smoothing  process  systematically  suppresses  the  eigenvectors 
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corresponding  to  eigenvalues  that  are  approximately  equal  to  the  current  value  of  the 
smoothing  parameter;  this  can  be  seen  by  observing  the  change  of  the  eigenvalues  printed 
out  by  the  code  at  each  value  of  the  smoothing  parameter. 

Twomey’s  criterion  for  the  degree  of  independence  of  the  measurements  can  be 
exploited  to  yield  an  indication  of  the  resolution  attainable  in  a  given  series  of 
measurements  (i.e.,  the  number  of  size  bins  which  can  be  accurately  deconvolved  in  a 
solution  histogram).  In  scattering  experiment  simulations  performed  at  AEDC,  15-,  30-, 
45-,  13  5-,  150-,  and  90-deg  angles  have  been  used;  both  polarization  states  have  been 
used,  and  the  particles  have  been  assumed  to  be  A1203,  illuminated  by  0.5145-p  light. 
For  these  conditions,  the  eigenvalue  criterion  indicated  that  good  deconvolutions  could 
be  obtained  for  a  maximum  of  six  bins,  and  the  simulated  deconvolutions,  even  including 
input  “noise”  obeyed  this  rule  of  thumb.  For  eight  bins,  non-physical  behavior  was 
manifested  by  the  tendency  of  the  solution  to  change  rapidly  with  smoothing  parameter 
from  a  form  having  negative  mole  fraction  values  to  an  asymptotic  “monotonic 
stairstep,”  declining  form.  Whenever  physical  information  is  present  in  the  solution,  there 
is  a  wide  range  of  smoothing  parameters  in  which  the  solution  exhibits  a  slow  transition 
between  the  extremes  mentioned  above. 

A  family  of  codes  denoted  as  SIZS1M  was  previously  constructed  on  the  basis  of  the 
constrained  linear  inversion  technique.  In  all  these  codes,  error  propagation  was 
considered  only  for  the  computation  of  residual  errors.  In  one  of  these  codes, 
experimental  limitations  prevented  the  use  of  the  90-deg  reference  signal.  In  this  code, 
the  reference  signal  was  taken  to  be  a  composite  of  the  unpolarized  signals  from  all  of 
the  angular  channels.  Results  of  this  analysis  are  reported  in  Ref.  A-4. 

A2.3  MODIFICATIONS  TO  INCLUDE  DECRAT  ERRORS 

In  constructing  code  S1ZFRED,  use  has  been  made  of  three  statistical  modifications 
of  the  basic  constrained  linear  inversion  technique.  These  modifications  permit  SIZFRED 
to  weight  the  DECRATS  inversely  proportional  to  the  level  of  imprecision  associated 
with  each  DECRAT.  In  addition,  upper  bound  output  imprecision  bands  associated  with 
each  size  bin  are  computed  from  a  statistical  error  propagation  analysis.  Finally, 
goodness-of-fit  parameters  have  been  devised  to  minimize  the  subjectivity  associated  with 
determination  of  the  optimum  smoothing  level  of  the  deconvolutions. 

The  first  modification  was  reported,  in  the  context  of  constrained  linear  inversion, 
in  Ref.  A- 10,  although  similar  techniques  are  discussed  in  standard  statistical  references 
such  as  Ref.  A-ll.  This  modification  is  implemented  in  SIZFRED  by  redefining  the 
symmetric  kernel  matrix  as  shown  below: 

Mil  -  l  Oki  %/(8Gl)2  (A'10) 
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In  Eq.  (A-IO)  the  set  of  quantities  |6Gkj denotes  the  input  imprecision  values  associated 
with  each  DECRAT,  and  the  input  errors  are  assumed  to  be  uncorrelated  (i.e.,  it  is 
assumed  that  the  input  covariance  matrix  has  only  diagonal  terms,  and  these  are  the 
squares  of  the  input  imprecision  values).  The  imprecision-weighted  solution  is 

Fi  -  S  Tk1  Ojk  Gj/tfGj)2  (A-!l) 

The  second  modification  of  the  basic  inversion  technique  is  also  a  simple 
generalization  of  standard  statistical  procedures  presented  in  Ref.  A-l  1 ;  however,  the 
implications  of  this  generalization  are  subtle  and  far  reaching,  as  discussed  in  detail  in 
Ref.  A-l 2. 

If  the  input  errors  are  assumed  to  be  statistically  independent,  then  the  upper 
bounds  on  the  output  imprecision  bands  can  be  computed  for  each  size  bin  from  the 
square  roots  of  the  variances  written  below: 

S2f,  -  Z*tfVSGk)2°llMH1  <A12) 

(Note  that  this  equation  is  not  summed  over  “i”.)  Upper  and  lower  bound  estimates  to 
the  solution  histogram  are  obtained  by  adding  the  square  roots  of  the  above  variances  to 
—  or  subtracting  these  quantities  from  —  the  solution  histogram  in  each  size  bin. 

Reference  A-l  2  shows  that,  because  of  the  use  of  constraints  in  the  deconvolution 
procedure,  the  extent  to  which  these  relations  overestimate  the  true  output  imprecision 
bands  cannot  be  assessed.  Therefore,  a  conservative  estimate  on  the  confidence  limits  is 
provided  with  which  the  histogram  solutions  should  be  associated.  For  example,  if  the 
DECRAT  imprecision  values  specified  as  inputs  to  SIZFRED  correspond  to  1 -sigma  values 
(68-percent  confidence  level),  the  code  will  yield  output  imprecision  bands  with  the 
following  nature:  The  probability  exceeds  0.68  that  solutions  which  have  random  input 
error  falling  within  the  stated  imprecision  bands  will  lie  within  the  computed  output 
imprecision  bands. 

In  reality,  there  is  little  certainty  that  an  input  Gaussian  error  distribution  will 
propagate  as  an  output  Gaussian  error  distribution.  Hence,  the  precise  meaning  of  the 
output  imprecision  bands  can  be  determined  only  by  performing  a  statistically  significant 
number  of  simulated  deconvolutions  of  computer  generated  scattering  data  with  various 
levels  of  Gaussian  distributed  error,  the  standard  distribution  of  which  corresponds  to  the 
imprecision  value  stated  for  each  DECRAT  and  the  mean  of  which  corresponds  to  each 
DECRAT.  Ideally,  at  least  20  such  computer  runs  should  be  made  for  each  of  the  three 
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types  of  size  distribution  functions  considered  here.  However,  such  a  comprehensive 
study  is  beyond  the  scope  of  this  report. 

The  third  statistical  modification  to  the  basic  deconvolution  technique  is  an  attempt 
to  determine  goodness-of-fit  parameters  that  will  reduce  the  subjectivity  associated  with 
determining  the  optimum  value  of  the  smoothing  parameter.  Using  strictly  heuristic 
arguments,  goodness-of-fit  parameters  have  been  derived  by  considering  the  propagation 
of  the  residual  error  and  including  the  magnitude  of  the  smoothing  parameter,  relative  to 
the  eigenvalues  of  the  smoothed,  symmetric  kernel  matrix.  Since  it  is  not  known  a  priori 
how  much  to  weight  each  eigenvalue,  the  goodness-of-fit  parameters  are  defined  as  shown 
below  for  two  different  ways  of  incorporating  the  eigenvalues: 

Ajt,B  =  |  RT/Qki^klRHlf(y)  (A- 13a) 

fA(y)  -  in 


fB  W  =  [n?=i la  ()A)]1/M  (a- 13c) 

Here  M  is  the  number  of  size  bins  in  the  histogram.  Equations  (A- 13  b)  and  (A- 13c), 
which  introduce  weighted  geometric  means  of  the  smoothing  parameter /eigen values  ratios, 
are  used  in  the  goodness-of-fit  parameters  respectively  denoted  GOOF  A  and  GOOF  B. 

A3.0  DETAILS  OF  CODE  SIZFRED 

A  flow  chart  of  code  SIZFRED  is  shown  in  Fig.  A-2.  The  program  requires  that 
precomputed  Mie  intensity  functions  be  read  in  at  equal  subintervals  of  size  parameter 
for  each  of  the  N+l  angles  at  which  experimental  data  are  to  be  entered.  These  functions 
must  have  the  form  shown  In  Eq.  (A-2).  The  code  assumes  that  the  last  angle  is  the 
reference  angle  (presumably  90  degrees).  The  tape  files  (which  are  generated  by  using 
code  MIESTEP)  must  fit  the  following  format:  The  first  variable  is  the  scattering  angle, 
the  second  is  the  size  parameter,  the  third  is  the  perpendicular  Mie  function,  and  the 
fourth  is  the  parallel  Mie  function.  It  is  assumed  that  the  Mie  functions  have  been 
integrated  over  the  experimental  range  of  subtended  scattering  angle.  If  the  trapezoidal 
rule  is  used  for  the  angular  integration,  numerical  instabilities  prevent  the  use  of 
scattering  kernels  with  subtended  scattering  angles  greater  than  about  one  degree.  Angular 
integration  by  means  of  Gaussian  quadrature,  however,  allows  larger  subtended  angles  to 
be  used,  and  it  is  recommended  that  the  kernels  be  computed  using  16-point  Gaussian 


(n.V,) 


1/M 


(A-I3b) 
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quadrature.  No  attempt  has  been  made  to  determine  the  maximum  permissible  subtended 
angles,  but,  clearly,  the  greatest  accuracy  is  obtained  by  limiting  the  subtended  angles  to 
the  smallest  values  consistent  with  adequate  counting  statistics  and  signal-to-noise  ratio. 

Table  A-l  shows  the  nomenclature  for  the  variables  that  the  code  user  requires,  and 
Table  A-2  shows  the  order  in  which  inputs  to  the  code  must  be  made.  Subroutine 
MIEFfLE-a  binary  search  routine -extracts  the  appropriate  Mie  functions  from  the  tape 
file  and  sends  them  to  the  numerical  integration  subroutines,  which  then  return  values 
from  which  the  elements  of  the  unsymmetric  scattering  kernel  matrix  are  constructed. 
MIEFILE  includes  a  variety  of  error  tests;  Table  A-3  defines  the  resulting  error  messages 
that  the  code  will  print  if  various  conditions  are  not  met. 

Once  the  unsymmetric  kernels  have  been  constructed,  the  code  symmetrizes  these 
matrices  and,  if  the  input  imprecision  values  have  been  stated,  weights  them  according  to 
the  DECRATS  having  the  least  imprecision.  The  symmetric  scattering  kernel  matrix  is 
then  added  to  the  product  of  the  smoothing  parameter  and  an  internally  generated 
smoothing  matrix,  and  the  resultant  smoothed  kernel  matrix  is  inverted  with  use  of 
Cholesky’s  method.  In  addition,  eigenvalues  of  the  smoothed  kernels  are  computed.  The 
matrix  inversion  routine  includes  a  test  to  determine  at  which  step  in  Cholesky’s 
factorization  process  significant  digits  may  have  been  lost.  This  test  is  necessary  because 
Cholesky’s  method  requires  that  the  matrix  being  inverted  be  positive  definite,  and  loss 
of  significant  digits  can  cause  the  matrix  to  lose  its  positive  definite  nature.  To  use  this 
error  check  correctly,  one  must  determine  the  error  criterion,  EPS,  by  trial  and  error. 
Generally,  a  value  on  the  order  of  0.000001  works  fairly  well, 

A  much  more  definitive  error  test  is  provided  by  multiplying  the  original  and 
inverted  matrices  and  by  then  testing  pairs  of  elements  along  the  diagonal.  When  the  code 
indicates  inaccuracies  in  the  nth  significant  digit,  off-diagonal  elements  of  the  product 
matrix  (which  should  be  the  identity  matrix)  have  magnitudes  on  the  order  of  the  n+lth 
significant  digit.  Thus,  the  inversion  error  criterion  is  conservative. 

After  matrix  inversion,  the  code  computes  the  solution  histogram,  the  residual  error 
vector,  the  norm  of  the  residuals,  and  the  relative  norm  of  the  residuals  (the  residuals’ 
norm  is  divided  by  the  DECRATS’  norm).  If  the  input  imprecision  values  have  been 
stated,  then  the  code  computes  output  imprecision  bands  and  prints  the  upper  and  lower 
bounds,  as  well  as  the  computed  histogram,  for  each  diameter  bin.  As  an  additional  check 
on  the  occurrence  of  computational  error,  the  code  also  computes  the  solution 
corresponding  to  zero  smoothing  level,  even  if  the  starting  value  of  the  smoothing 
parameter  is  nonzero.  Finally,  the  two  goodness-of-fit  parameters,  GOOF  A  and  GOOF  B, 
are  computed  and  printed.  This  process  is  repeated  until  value  of  the  smoothing 
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parameter  equals  or  exceeds  its  stop  value,  EPST.  At  this  point  the  computer  generates  a 
histogram  plot.  The  last  card  of  the  input  deck  specifies  the  plot  parameters.  The  first 
variable  on  this  card  is  the  minimum  diameter  that  will  be  plotted  and  the  second  is  a 
scale  factor.  For  six-bin  histograms  this  factor  should  be  0.25. 

A4.0  ILLUSTRATIVE  SIZING  CALCULATIONS 

Simulated  PSDF  deconvolutions  have  been  performed  on  the  three  types  of 
distribution  functions  shown  in  Figs.  A-3  through  A-5.  The  polarized  angular  scattering 
distributions  corresponding  to  the  size  distribution  of  Fig.  A-3  are  shown  in  Fig.  A-6.  The 
corresponding  DECRATS  for  the  previously  given  scattering  angles  are  listed  in  Tables 
A-4  through  A-6  for  the  three  size  distributions  shown  in  Figs.  A-3  through  A-5. 

When  a  deconvolution  of  the  data  in  Table  A-4  is  attempted  without  specifying  the 
(arbitrarily  chosen,  but  presumably  typical)  set  of  imprecision  values  stated  in  Table  A-4, 
no  minimum  is  found  in  either  GOOF  parameter.  Plotting  these  parameters  against  the 
smoothing  parameter.  7,  however,  reveals  that  a  value  in  the  range  from  0.004  to  0.006 
should  be  nearly  optimum;  the  solution  changes  negligibly  within  this  range.  The  resulting 
deconvolution  is  shown  in  Fig.  A-7  and  can  be  compared  with  the  original  size 
distribution  function  in  Fig.  A-3  by  dividing  the  computed  size  histogram  mole  fractions 
by  the  histogram  bin  size. 

The  result  of  deconvolution  involving  the  stated  imprecision  values  (of  Table  A-4)  is 
shown  in  Fig.  A-8.  The  histogram  corresponds  to  a  minimum  in  the  GOOF  B  parameter, 
which  occurred  when  the  smoothing  parameter  just  exceeded  the  smallest  eigenvalue  of 
the  smoothed  kernel  matrix.  Because  of  the  dominance  of  the  scattering  functions  by 
larger  particles,  the  imprecision  of  the  smallest  size  bin  is  rather  large.  A  dramatic 
improvement  occurs,  however,  when  the  magnitude  of  all  the  input  imprecision  bands  is 
halved,  as  shown  in  Fig.  A-9,  Similar  but  even  more  dramatic  improvement  in  output 
imprecision  is  shown  in  the  deconvolution  of  the  simulated  bimodal  distribution  (Figs. 
A-4,  A- 10,  and  A-ll).  The  DECRATS  corresponding  to  this  case  were  obtained  by 
averaging  the  scattering  functions  obtained  from  two  narrow  Deirmendjian  distribution 
functions.  A  previous  attempt  to  resolve  bimodal  distributions  with  closely  spaced  peaks 
was  unsuccessful.  Minima  in  both  GOOF  parameters  were  found  in  the  successful  bimodal 
deconvolutions.  The  solution  histograms  differed  very  little  for  the  two  minima. 

Figures  A-5  and  A-12  are  included  to  illustrate  the  ability  of  the  deconvolution 
procedure  to  treat  a  very  narrow  distribution  function.  It  was  feared  that  the  technique 
would  attempt  to  “fill  up”  empty  bins,  but  Fig.  A-12  shows  that  this  did  not  happen. 
For  the  smoothing  parameter  increment  used,  it  was  not  possible  to  find  a  solution  in 
which  there  were  no  negative  mole  fractions  (though  the  magnitude  of  the  negative  mole 


55 


AEDC-TR-80-3 


fractions  is  quite  small).  Although  Fig.  A-12  does  correspond  to  an  optimum  smoothing 
level,  it  is  likely  that  a  smaller  increment  of  smoothing  parameter  could  have  located 
histograms  near  the  optimum  smoothing  level  for  which  only  positive  mole  fraction 
values  occurred.  Nevertheless,  Fig.  A-12  shows  that  SIZFRED  will  not  fill  up  empty  size 
bins  unless  the  smoothing  parameter  is  allowed  to  become  so  large  that  clearly  unphysical 
results  occur. 

As  previously  stated,  it  is  beyond  the  scope  of  this  work  to  perform  a  detailed  test 
of  the  propagation  of  error  distribution  shapes  through  the  deconvolution  procedure. 
However,  a  few  test  cases  have  been  investigated  that  had  simulated  error  included  in  the 
DECRATS.  Figure  A-13  shows  the  deconvolution  resulting  from  simulating  error 
occurrence  only  in  DECRAT  1,  as  shown  in  Table  A-7.  A  more  stringent  test  is  afforded 
by  simulating  Gaussian  distributed  errors  as  shown  in  Table  A-8.  For  this  set  of  errors 
(corresponding  to  the  imprecision  values  in  Table  A-4),  the  eigenvalue  criterion  shows 
that  a  six-bin  deconvolution  cannot  be  obtained;  this  is  confirmed  by  the  behavior  of  the 
solution  histogram  as  7  is  increased.  The  same  criterion  indicates  that  a  five-bin  histogram 
can  be  obtained.  The  result  corresponding  to  a  minimum  in  GOOF  A  is  shown  in  Fig. 
A-14.  Better  results  are  obtained  by  increasing  the  smoothing  parameter  until  the  two 
GOOF  parameters  agree  to  within  about  20  percent,  as  shown  in  Fig.  A-l  5.  However,  this 
result  must  be  regarded  as  fortuitous  until  more  extensive  tests  with  other  distribution 
functions  and  other  error  distributions  can  be  conducted. 

Decreasing  the  input  imprecision  by  a  factor  of  two  and  recomputing  a  Gaussian 
error  distribution  (Table  A-9)  result  in  a  six-bin  deconvolution  that  is  marginally  allowed 
by  the  eigenvalue  criterion  (Fig.  A-l 6).  Again,  the  histogram  corresponding  to  near 
equality  of  the  GOOF  parameters  is  superior  to  that  corresponding  to  a  minimum  in 
either  parameter. 


A5.0  CONCLUDING  REMARKS 

The  primary  objectives  of  this  study  were  to  incorporate  error  propagation  analysis 
into  AEDC’s  size  deconvolution  codes  and  to  develop  a  means  of  identifying  the 
optimum  smoothing  level  with  as  little  subjectivity  as  possible.  Code  SIZFRED  has  been 
developed  to  accomplish  these  objectives.  Although  some  judgement  must  still  be 
exercised  in  determining  the  optimum  smoothing  level  and  in  determining  whether  the 
results  are  physically  meaningful,  progress  toward  minimizing  the  subjectivity  has  been 
achieved  by  the  introduction  of  two  goodness-of-fit  parameters.  Extensive  testing  of  the 
code  must  yet  be  undertaken  to  determine  the  extent  of  transformation  of  a  Gaussian 
error  distribution  by  error  propagation  and  to  determine  how  such  error  propagation 
transformations  affect  the  characteristics  and  usefulness  of  the  existing  error  criteria. 
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The  effects  of  other  uncertainties,  such  as  particle  nonsphericity  and  complex 
refractive  index  uncertainty,  have  not  been  addressed.  Ideally,  the  deconvolution 
procedure  should  be  tested  further  by  simulating  polydisperse  multispecies  scattering, 
perhaps  using  approximate  nonspherical  scattering  functions.  The  code  should  also  be 
tested  in  laboratory  scattering  experiments  with  polydispersions  of  scatterers  with  known 
properties,  such  as  commercially  available  latex  spheres. 

As  configured  presently,  code  S1ZFRED  can  accept  inputs  from  21  scattering  angles 
(in  both  polarization  states)  and  will  produce  deconvolutions  with  up  to  30  size  bins 
(depending  on  angular  detector  locations,  error  levels,  etc.).  Whether  such  experimental 
complexity  is  justified,  in  terms  of  attainable  size  resolution,  has  not  been  addressed. 
Parametric  studies  should  be  carried  out  with  various  size  distributions  and  varying  error 
levels  to  attempt  to  optimize  the  number  of  detector  channels  needed  and  their  angular 
locations. 

These  preliminary  results  show  that  the  deconvolution  procedure  is  stable  for  errors 
as  large  as  ±3  to  5  percent  in  each  DECRAT  and  at  the  expense  of  loss  of  histogram 
resolution  (i.e.,  number  of  bins).  Moreover,  if  the  input  errors  are  restricted  to  ±1  to  3 
percent,  then  these  studies  also  suggest  that  deconvolution  can  be  achieved  without 
significant  loss  of  size  resolution  and  with  substantial  decrease  in  output  imprecision  in 
small  particle  size  mole  fractions.  To  the  extent  that  it  is  practical,  experimental 
applications  should  be  directed  toward  achieving  these  detection  imprecision  levels. 
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Figure  A-2.  How  chart  of  code  SIZFRED. 
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Figure  A-4.  Bimodal  size  distribution  function  for  t8St  of 
deconvolution  procedure. 
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Figure  A-5.  Center-peaked  size  distribution  function  for  test  of 
deconvolution  procedure. 
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Figure  A-7.  Deconvolution  histogram  obtained  from  DECRATS  in  Table  A-4 
without  using  stated  input  imprecision  values  {corresponds 
to  standard  size  distribution  function). 
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Figure  A-8.  Deconvolution  histogram  obtained  from  DECRATS  in  Table  A-4 
using  stated  input  imprecision  values  (corresponds 
to  standard  size  distribution  function). 
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Figure  A-9.  Deconvolution  histogram  obtained  from  DECRATS  in  Table  A-4 
using  %  stated  input  imprecision  values  (corresponds 
to  standard  size  distribution  function). 
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Figure  A-10.  Deconvolution  histogram  obtained  from  DECRATS  in  Table  A-5 
using  stated  Imprecision  values  (corresponds  to 
bimodal  size  distribution  function]. 
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Figure  A-1 1 .  Deconvolution  histogram  obtained  from  DECRATS  in  Table  A-5 
using  %  stated  input  imprecision  values  (corresponds  to 
bimodal  size  distribution  function). 
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Figure  A-1 2.  Deconvolution  histogram  obtained  from  DECRATS  in  Table  A-6 
using  stated  input  imprecision  values  (corresponds  to 
center-peaked  distribution). 


Figure  A-1 3.  Deconvolution  histogram  obtained  from  DECRATS  in  Table  A- 7 
using  stated  input  imprecision  values  (corresponds  to 
standard  size  distribution  function  with  8.6-percent 
error  inserted  into  DECRAT  1  only). 
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Figure  A- 14.  Five-bin  deconvolution  histogram  obtained  from 

DECRATS  in  Table  A-8  using  optimum  solution  based  on  minimum 
in  GOOF  A  {corresponds  to  standard  size  distribution  function 
with  Gaussian  error  distribution  inserted  into  all  DECRATS) 


Figure  A-1 5.  Five-bin  deconvolution  histogram  obtained  from 

DECRATS  in  Table  A-8  using  optimum  solution  based  on  equivalence 
of  GOOF  A  and  GOOF  B  to  20-percent  {corresponds  to  standard 
size  distribution  function  with  Gaussian  error  distribution 
inserted  into  all  DECRATS) 


Figure  A-1 6.  Six-bin  deconvolution  histogram  obtained  from  DEGRATS  in 
Table  A-8  using  optimum  solution  based  on  equivalence  of 
GOOF  A  and  GOOF  B  to  20-percent  (corresponds  to 
standard  size  distribution  function  with  Gaussian  error 
distribution  inserted  into  all  DECRATS). 
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Table  A-1 .  Nomenclature  Used  in  Code  SIZFRED 


THETA  (I);  1=1, M 
THETAR  =  Theta  (n+1) 
DMAX 

DMTN 

HSIZ 


HD  I  AM 

SIZMAX 

SIZMIN 

DELBIG 

_  SIZHAX-SIZMIN 
L  “  HSIZ 

M 

DELSIZ 


NDIM 


DELSIZ 

HSIZ 


Scattering  angle,  deg 
Reference  angle,  deg 

Maximum  particle  diameter,  y,  considered 
for  histogram 

Minimum  particle  diameter,  y,  considered 
for  histogram 

Size  parameter  increment  at  which  Mie 
functions  have  been  recorded  on  tape 
called  by  code  SIZFRED 

Diameter  increment,  j,  corresponding 
to  HSIZ 

Size  parameter  corresponding  to  DMAX 
Size  parameter  corresponding  to  DMIN 
SIZMAX  -  SIZMIN 

Number  of  size  parameter  increments 
between  SIZMIN  and  SIZMAX 

Number  of  bins  in  histogram 

Size  parameter  range  corresponding  to 
one  bin 

Number  of  size  parameter  increments 
in  one  bin 


DEL 


DMAX-DMIN 

M 


Diameter  range,  y,  corresponding 
to  one  bin 


D(J);  J=1 ,M+L 


S(K) ;  K=L,L 


A(I,J);  I-1.2H; 
J=1  ,M 


AREF 

Q(I,J)  = 

1=1, 2N; 

AM(I,J); 

J=l,M 


A(I,J). 
AREF  ’ 

J=1  ,M 

i=l. M; 


Diameters,  y,  spaced  at  intervals 
of  DEL  from  DMIN  to  DMAX 

Size  parameters  spaced  at  intervals  of 
HSIZ  from  SIZMIN  to  SIZMAX 

Mie  functions  averaged  over  DELSIZ, 
returned  by  subroutine  INTEGS 

Reference  Mie  functions  averaged  over 
DELBIG,  returned  by  subroutine  INTEGR 

Unsymmetric  scattering  kernel  matrix 


Symmetric  scattering  kernel  matrix 
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I-i.M; 

J=1,M 

R(i,J); 

I-1.H; 

J-1,M 

AR(I,J) 

;BR(1,J) 

G(K);K=1,2N 

DG(K);K=1,2N 

BPPS(K);  K=l,2N 

ERR 

ERREL 

GAMO 

DELGAM 

EPST 

EPS 

IER 


P(I);  1=1, M 


FH(I) ;  1=1, M 

AM 

CN 

V(I)S  I-1.M 
DELF(I);  1=1, M 

GOOF  A 


GOOF  B 


Table  A-1.  Concluded 

Smooching  matrix 

Smoothed,  symmetric  kernel  matrix 
(Mote:  This  matrix  Is  redefined  In 
calculation  of  goodness-of-fit  parameter.) 

Dumny  matrices  used  to  prevent  destruction 
of  R(I,J)  in  inversion  and  eigenvalue 
calculations.  Inverse  matrix  Is  returned 
In  place  of  AR;  diagonalized  matrix  Is 
returned  In  place  of  BR. 

Vector  of  input  deconvolution  ratios 
(DECRATS) 

Vector  of  input  imprecision  estimates 

Vector  of  residual  errors 

Norm  of  EPPS(K)  vector 

Ratio  of  ERR  to  norm  of  G(K)  vector 

Initial  value  of  smoothing  parameter 

Smoothing  parameter  increment 

Final  value  of  smoothing  parameter 

Input  constant  used  in  accuracy  check 
in  Cholesky’s  matrix  inversion  routine 

Error  parameter  returned  if  loss  of 
significant  digits  occurs  in  matrix 
inversion 

Solution  histogram  (mole  fraction  in 
each  bin) 

Unsmoothed  solution  histogram 

Normalization  constant  for  F(I) 

Normalization  constant  for  FH(I) 

Solution  imprecision  bands 

Temper ary  vector  used,  separately,  in 
computation  of  FH(I)  and  in  determina¬ 
tion  of  goodness-of-fit  parameters 

Goodness-of-fit  parameter  based  on 
propagation  of  residual  error  and  on 
the  logarithm  of  the  ratio  of  the 
smoothing  parameter  to  the  geometric 
mean  of  all  kernel  eigenvalues 

Goodness-of-fit  parameter  based  on 
propagation  of  residual  error  and  on 
geometric  mean  of  the  logarithms  of 
the  ratios  of  the  smoothing  parameter 
to  each  kernel  eigenvalue 
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Table  A-2.  SIZFRED  Inputs 


Card  Number 

Parameter 

Columns 

1. 

Name  for  plot  cover  sheet 

1  through  16 

2. 

Title  for  plot  and  prlnt-out 

1  through  72 

3.  a) 

b) 

c) 
d> 

Initial  smoothing  parameter  GAMO 

Smoothing  parameter  Increment  DELGRAM 

Final  smoothing  parameter  EPST 

Tolerance  parameter  used  in  determining 
step  at  which  matrix  inversion  lost 
accuracy  —  EPS 

1  through  LO 

1 1  through  20 

2 1  through  30 

3 1  through  40 

4. a) 

b) 

c) 

Initial  size  parameter  —  SIZMIN 

Light  wavelength  in  microns  —  XLAM 

Size  parameter  increment  —  HSIZ 

1  through  10 

11  through  20 

21  through  30 

5. a) 

Number  of  size  parameter  increments 
between  minimum  and  maximum  size 
parameters  —  L 

1  through  4 
Right  Justified 

b) 

Number  of  scattering  angles  (excluding 
reference  angle)  at  which  data  will  be 
entered  —  N 

5  through  8 
Right  Justified 

c) 

Number  of  size  parameter  increments  in 
each  size  bin  —  NDIM 

9  through  12 
Right  Justified 

6. a) 

Deconvolution  ratios  G(I);  1-1, 2N 

1  through  10 

b) 

Imprecision  values  DG(I);  1=1 P2N 

There  are  2N  cards  of  this  type. 

11  through  20 

7. 

Scattering  angles  (in  degrees) 

THETA(I);  I-l.N+l 

There  are  NH  cards  of  this  type. 

The  last  card  specifies  the 
reference  angle 

1  through  10 

8. a) 

Minimum  diameter  plotted  on  histogram 
(microns) 

1  through  10 

b> 

Scale  factor  for  histogram 

11  through  20 
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Table  A-3.  Error  Messages  Generated  by  Subroutine  Mlefile 

Flag 

0  Everything  is  OK 

1  Either  the  angle  THETA  was  not  found  or  the  size 
parameters  sought  were  larger  than  the  largest 
size  parameter  in  the  files 

2  The  first  size  parameter  was  smaller  than  the 
smallest  in  the  files 

3  The  correct  angle  and  range  of  size  parameters 
were  found  but  the  first  size  parameter  did  not 
coincide  with  any  of  the  size  parameters  in  the 
file 

4  The  correct  angle,  range  of  size  parameters,  and 
the  first  size  parameter  were  found,  hut  some 
other  size  parameter  sought  was  not  found  in  the 
file 

5  Same  as  above,  but  the  search  continued  into 
subsequent  files 

6  Same  as  above,  but  the  search  ran  into  the  end 
of  the  tape 


Table  A-4.  Deconvolution  Ratios  and  Input  Imprecision  Values  Corresponding 
to  the  Standard  Size  Distribution  Function  (Shown  in  Fig.  A-3) 


DECRATS 

Standard  Deviations 

6.894 

0.6 

3.203 

0.3 

2.388 

0.2 

0.3657 

0.04 

0.4197 

0.04 

7.093 

0.3 

3.082 

0.2 

2.3805 

0.2 

0.5769 

0.06 

0.451 

0.1 

Relative  Standard  Deviation  =  0.0722 
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Table  A-5. 


Table  A-6. 


Deconvolution  Ratios  and  Input  Imprecision  Values  Corresponding 
to  the  Bimodal  Size  Distribution  Function  (Shown  in  Fig.  A-4) 

DECKATS  Standard  Deviations 


9.157 

0.9 

3.566 

0.4 

2.672 

0.3 

0.3305 

0.03 

0.4224 

0.04 

9.723 

0.9 

3.578 

0.4 

2.721 

0.3 

0.5192 

0.05 

0.3325 

0.03 

Relative  Standard  Deviation  =  0.0985 


Deconvolution  Ratios  and  Input  Imprecision  Values  Corresponding 
the  Center-Peaked  Size  Distribution  Function  (Shown  in  Fig.  A-5} 


DECRATS 

Standard  De1 

8.017 

0.6 

2.771 

0.3 

2.227 

0.2 

0.3744 

0.04 

0.4486 

0.03 

8.157 

0.6 

2.433 

0.1 

1.994 

0.2 

0.504 

0.05 

0.406 

0.04 

Relative  Standard  Deviation  =  0.0767 


74 


AEDC -TR-80-3 


Table  A- 7.  Deconvolution  Ratios  and  Input  Imprecision  Values  Corresponding 
to  the  Standard  Size  Distribution  Function  but  Having 
8.6-percent  Error  in  DECRAT  1 

DECRATS  Standard  Deviations 


6.3 

0.6 

3.203 

0.3 

2.388 

0.2 

0.3697 

0.04 

0.4197 

0.04 

7.093 

0.3- 

3.082 

0.2 

2.3805 

0.2 

0.5769 

0.06 

0.451 

0.1 

Relative  Standard  Deviation  =  0.0745 


Table  A-8.  Original  and  Gaussian  Distributed  Deconvolution  Ratios 
(Corresponding  to  the  Standard  Deviations  in  Table  A-4) 

DECRATS  Gaussian  Distributed  Vectors 


6.894 

6.715 

3.203 

2.917 

2.388 

2.360 

0.3657 

0.3629 

0.4197 

0.4377 

7.093 

6.811 

3.082 

3.289 

2.3805 

2.207 

0.5769 

0.5281 

0.451 

0.3633 

Relative  Error  Level  -  0.0463 
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Table  A- 9.  Original  and  Gaussian  Distributed  Deconvolution  Ratios 
[Corresponding  to  One-Half  the  Standard  Deviations  in 
Table  A-4) 


DECRATS 

6.894 

3.203 

2.388 

0.3657 

0.4197 

7.093 

3.082 

2.3805 

0.5769 

0.451 


Gaussian  Distributed  Vectors 
6.805 
3.060 
2.374 
0.3643 
0.4287 
6.952 
3.185 
2.294 
0.5526 
0.4072 


Relative  Error  Level  =  0.0231 
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APPENDIX  B 

MEASUREMENT  OF  PARTICLE  SIZES  AND 
COMPLEX  REFRACTIVE  INDICES  BY  USE  OF 
MIE  SCATTERING  TECHNIQUES 


(Presented  at  Smoke  Symposium  II,  Harry  Diamond  Laboratory, 
Adelplii,  Maryland,  April  25  -26,  1978). 

B.  P.  Curry,  D.  P.  Weaver,  and  J.  W.  L.  Lewis. 

ARO,  Ine. 
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Recently,  an  experiment  was  carried  out  to  determine  the  size  distribution,  number 
density,  and  chemical  identity  of  particles  present  in  the  flow  of  the  Arnold  Engineering 
Development  Center  (AEDC>von  Karmun  Gas  Dynamics  Facility  (VKF)  Tunnel  F,  a 
hypervelocity,  arc-drive,  intermittent-flow  wind  tunnel.  This  appendix  describes  the 
results  of  that  experiment  and  the  techniques  used  in  the  analysis.  In  addition,  a 
procedure  is  described  for  obtaining  the  complex  indices  of  refraction  of  particles  at 
wavelengths  at  which  these  indices  are  uncertain. 

Copper  particles  and/or  tungsten  particles  were  suspected  to  be  present  in  the  flow, 
since  these  materials  comprise  the  electrodes  between  which  the  arc  is  struck.  Also,  iron 
particles  might  have  been  present  because  steel  baffles  and  burst  diaphragms  were  used. 
In  addition,  the  expansion  of  the  flow  through  the  hypersonic  nozzle  was  expected  to 
produce  condensation  of  both  water  vapor  and  nitrogen  present  in  the  ambient  air. 
Although  all  of  the  particles  were  potential  sources  of  the  observed  light  scattering, 
subsequent  analysis  showed  that  the  data  could  be  explained  in  only  two  ways:  (I)  of 
single  species  present,  only  copper  could  explain  the  scattering  data  and  (2)  of  dual 
species  combinations,  copper  and  tungsten  with  no  more  than  2 1 -percent  mole  fraction 
(44-percent  mass  fraction)  of  tungsten  present  could  explain  the  analysis  of  the  chosen 
data  point  -  which  was  in  the  middle  of  a  1 50-msec  run.  Later,  chemical  analysis  of  the 
mass  loss  from  the  electrodes  during  the  entire  run  showed  the  overall  mass  fraction  to 
be  50-percent  tungsten.  Analysis  of  other  data  points  during  the  run  is  currently  in 
progress.  The  tungsten/copper  ratio  is  expected  to  change  progressively  during  the  run. 

Figure  8-1  shows  the  experimental  configuration.  The  particle  flow  is  from  left  to 
right  at  about  105  cm/sec  speed.  A  pulsed  argon  ion  laser  (0.5145  wavelength) 
polarized  at  45  deg  to  the  scattering  plane  illuminated  the  particles.  A  detector  was 
located  opposite  the  laser  to  monitor  transmission,  and  two  spectrometers  were  used  to 
monitor  scattered  light.  (During  the  course  of  the  experiment  it  became  obvious  that 
spectrometers  were  unnecessary  and  that  ordinary  phototubes  would  work  just  as  well.) 
A  chopper  was  interposed  in  each  scattering  channel,  the  plan  of  which  is  shown  in  Fig. 
B-2.  Each  chopper  was  divided  into  groups  of  three  segments.  The  first  segment  in  each 
group  was  blank,  to  provide  a  marker  pulse  of  unpolarized  light.  The  second  and  third 
segments  of  each  group  consisted,  respectively,  of  vertically  and  horizontally  oriented 
polarizers.  The  chopper  speed  was  such  as  to  provide  10-msec  time  resolution.  Thus, 
about  10  data  points  could  be  obtained  during  the  “good  flow”*  portion  of  each  short 
run.  The  transmission  data  indicated  that  only  single  scattering  was  occurring  during  the 
good  flow  period,  though  multiple  scattering  clearly  occurred  prior  to  the  establishment 
of  the  flow  and  after  flow  breakdown,  as  indicated  by  the  nearly  complete  obscuration 


♦The  “good  How”  period  is  defined  by  “smooth  decline”  of  the  tunnel  pressure  versus  time. 
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of  the  transmitted  light.  Because  of  geometrical  constraints  imposed  by  the  tunnel 
observation  ports,  the  scattering  angles  used  were,  respectively,  146.5  deg  with  subtended 
angles  of  2.65  deg,  and  161.9  deg  with  subtended  angles  of  6.6  deg.  Subsequent 
computer  simulations  have  determined  that  better  accuracy  can  be  achieved  by  use  of 
smaller  subtended  angles,  and  future  experiments  will  be  planned  accordingly. 

Average  particle  size,  species  identification,  and  size  histograms  were  determined 
from  various  scattering  ratios  in  order  to  avoid  the  difficulties  associated  with  absolute 
intensity  measurements  involving  unknown  particle  number  density.  The  first  of  these 
ratios  is  the  degree  of  polarization  shown  (for  copper)  for  the  two  scattering  channels  in 
Figs.  B-3  and  B-4.  The  theoretical  curves  shown  have  been  numerically  integrated  over 
the  appropriate  subtended  angles.  The  experimental  imprecision  bunds  correspond  to  one 
standard  deviation  in  the  scattered  intensities.  Figures  B-3  and  B-4  eliminate  much  of  the 
non-uniqueness  in  the  determination  of  average  particle  size.  In  addition,  the  transmission 
data  can  be  combined  with  the  scattering  data  to  provide  "angular/exti  action”  ratios  - 
shown  in  Fig.  B-5  for  a  single  polarization  state  of  146.5  deg.  These  ratios  arc 
proportional  to  the  Mie  phase  functions  for  each  separate  polarization  state. 

These  ratios  are  obtained  experimentally  as  follows:  (l)  From  the  transmission  data, 
one  obtains  the  product  of  particle  number  density  times  extinction  cross  section  times 
length  over  which  scattering  and  absorption  occur,  (2)  From  the  angular  scattering  data, 
one  obtains  the  product  of  number  density  times  scattering  cross-section  times  length  of 
scattering  volume  (determined  entirely  by  the  collection  optics).  Thu  total  length  over 
which  extinction  occurs  could  only  be  estimated  to  lie  be L ween  the  nozzle  exit  diameter 
and  the  flow-field  diameter  obtained  from  method-of-characteristics  calculations.  The 
overlapping  experimental  imprecision  bands  shown  in  Fig.  B-5  were  thus  obtained.  For 
each  particle  species  there  are  four  different  plots  of  the  type  shown  in  Fig.  B-5. 
Comparison  of  these  results  with  the  results  oT  the  polarization  analysis  completely 
removes  the  non-uniqueness  in  the  average  size  and  shows  that  the  average  diameter  is 
0.8  ±  0.1  micron.  Using  the  calculated  extinction  cross  section  in  conjunction  with  the 
transmission  data  shows  the  particle  number  density  to  be  equal  to  2  x  104  cm*3  ±  10 
percent.  The  limiting  factor  in  average  size  determination  is  the  accuracy  with  which  the 
extinction  length  can  be  determined. 

The  above  procedure  was  repeated  for  the  other  particle  species,  and  all  species 
except  copper  and  tungsten  were  systematically  eliminated.  During  this  analysis,  it 
became  apparent  that  the  degree  of  polarization  is  insensitive  to  the  imaginary  part  of  the 
refractive  index  but  is  very  sensitive  to  the  real  part  of  the  index  of  refraction,  especially 
near  the  “rainbow  angle”  (which  is  not  very  sensitive  to  size).  In  contrast,  the  magnitude 
of  the  phase  function  plots  is  very  sensitive  to  the  imaginary  index,  but  relatively 
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insensitive  to  the  real  part.  Hence,  the  following  procedure  is  proposed  for  determining 
unknown  complex  refractive  indices  of  aerosols  in  the  infrared  spectral  region: 

1 .  With  use  of  the  visible  wavelengths  at  which  the  refractive  indices  are 
known  for  a  given  aerosol  species,  the  average  size  and  size  distribution 
histograms  are  determined  (by  the  method  to  be  described  shortly.) 

2.  The  real  refractive  index  is  iteratively  determined  by  measuring  degree  of 
polarization  in  the  infrared  at  angles  near  the  suspected  rainbow  angle  and 
by  then  adjusting  the  real  index  until  polydisperse  Mie  calculations  made 
with  the  previously  determined  particle  sizes  agree  with  the  experimental 
data. 

3.  Similarly,  the  imaginary  index  is  determined  by  adjusting  it  until  the 
polydisperse  Mie  phase  functions  agree  with  the  experimental  results. 

This  procedure,  while  admittedly  somewhat  non-unique,  should  yield  valuable 
information  if  it  is  judiciously  applied. 

Our  procedure  for  determining  particle  size  histograms  is  a  modification  of 
Twomey’s  constrained  linear  inversion  procedure.  Computer  simulations,  such  as  that  of 
Fig.  B-6,  have  shown  that  accurate  deconvolutions  can  be  obtained  for  a  limited  number 
of  sizeN  bins  by  deconvolving  ratios  of  angular  scattered  intensities  in  separate  polarization 
states  to  the  unpolarized  90-deg  scattered  signal.  In  Fig.  B-6,  the  stated  polydisperse 
distribution  has  been  integrated  over  the  size  bins  shown  -  the  results  being  the  dashed 
line  histogram.  The  solid  line  histogram  is  the  result  of  deconvolving  computer  simulated 
scattering  measurements.  The  accuracy  indicated  requires  subtended  angles  no  larger  than 
1  deg.  In  addition,  at  the  stated  small  refractive  indices*  and  the  indicated  forward 
scattering  angles,  no  more  than  three  bins  could  be  deconvolved  accurately.  Larger  indices 
and/or  larger  scattering  angles  permit  more  bin  deconvolutions. 

Because  of  the  unavailability  of  a  90-deg  reference  channel,  the  experimental 
deconvolutions  were  obtained  as  follows:  The  angular  scattering  intensities  for  each  angle 
and  each  polarization  state  were  ratioed  to  the  sum  of  the  data  taken  at  both  angles  and 
at  both  polarization  states.  These  ratios  were  then  used  as  inputs  into  the  constrained 
linear  inversion  program,  and  the  smoothing  parameter  was  adjusted  until  the  computed 
average  diameter  approximately  equaled  the  average  diameter  determined  previously. 
Computer  simulation  showed  three-bin  deconvolutions  obtained  in  this  fashion  to  be 
quite  accurate.  The  four-bin  deconvolutions  proved  to  be  marginal  by  Twomey’s  criterion 


•Deliberately  chosen  as  a  “worst  case." 
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regarding  the  size  of  the  smallest  eigenvalue  of  the  scattering  kernel  matrix  relative  to  the 
magnitude  of  the  square  error.  Nevertheless,  as  Fig.  B-7  shows,  the  experimental 
deconvolutions  are  reasonable  in  comparison  with  the  deconvolution  of  a 
computer-simulated  polydispersion  which  yields  about  the  same  average  diameter. 

The  limits  of  these  techniques  are  currently  being  established  at  AEDC  in  a  series  of 
laboratory  experiments  and  in  other  computer-simulation  studies.  All  indications  to  date 
are  that  these  techniques  should  prove  to  be  a  valuable  tool  in  establishing  aerosol  size 
and  optical  property  data. 
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Figure  B-1 .  Schematic  diagram  of  apparatus  placement. 
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b.  Partial  view  ol  chopper 


Figure  B-2.  Schematic  of  components  in  optical  path. 
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Figure  B-3.  Degree  of  polarization  at  146.5  deg. 
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Figure  B-4.  Degree  of  polarization  at  161 .9  deg. 
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Figure  B-5.  Typical  angular/extinction  ratio  plot. 
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Figure  B-6.  Deconvolution  of  a  simulated  polydispersion. 
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Figure  B-7.  Comparison  of  deconvolved  sizing  measurements  with 
simulated  polydispersion  deconvolution. 
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NOMENCLATURE 

(For  Main  Text  Only) 

i,  Perpendicular  scattering  ratio  for  ith  detector  channel 

i,  ^i)/$o  Parallel  scattering  ratio  for  ith  detector  channel 
Particle  diameter 

Normalized  particle  size  distribution  function,  PSDF 
Matrix  defined  in  Eq.  (19) 

Mie  angular  scattering  intensity  functions  for  particles  of 
.diameter  D  and  complex  refractive  index 

Matrix  defined  in  Eq.  (19) 

Propagation  wave  number  of  incident  radiation 

Number  of  subintervals  in  the  size  parameter  interval 

xm  in  ^  ^  xm  ax 

Line-of-sight  extinction  distance  in  Eq.  (3) 

Scattering  length  of  the  focal  volume  [i.e.,  projection  of  the 
focal  volume  (subtended  by  the  collection  optics  associated  with 
any  given  detector)  along  the  incident  laser  axis] 

Function  defined  by  Eq.  (24) 

Number  of  scattering  angles  at  which  measurements  are  made 

Number  density  of  particles  (i.e.,  number  of  particles  per  unit 
volume) 

Degree  of  polarization  of  scattered  radiation  [Eqs.  (6)  through 
(9)] 

Function  defined  by  Eq.  (24) 

Efficiency  factor  (dimensionless)  for  scattering,  absorption,  or 
extinction  by  particles  of  diameter  D  and  complex  refractive 
index  tj  =  rj,  -  iK 
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Angular/Extinction  Ratios  (ANGER  functions)  for  perpendicular 
and  parallel  polarization  states,  respectively 

Angular/Backscattcring  Ratios  (ANGBAR  functions)  for 
perpendicular  and  parallel  polarization  states,  respectively 

Ratio  of  transmitted  radiation  flux  to  incident  flux 

Matrix  defined  in  Eq.  (21) 

Size  parameter  for  particle  of  diameter  D  (dimensionless) 

Particle  chemical  species  mole  fraction 

Range  of  scattering  angles  subtended  by  the  collection  optics  of 
any  given  detector 

Range  of  azimuthal  angles  subtended  by  the  collection  optics  of 
any  given  detector 

Solid  angle  subtended  by  the  collection  optics  of  any  given 
detector 

Estimated  imprecision  in  a  measured  quantity 

Tolerance  parameter  used  in  multispecies  average  size 
determination 

Complex  refractive  index  of  particles 

Real  (or  dispersive)  refractive  index  component 

Polar  scattering  angle  of  spherical  coordinate  scattering  geometry 
shown  in  Fig.  5 

Imaginary  (or  absorptive)  refractive  index  component 

Wavelength  of  incident  radiation 

Cross  section  for  scattering,  absorption,  or  extinction 

Differential  cross  section  per  unit  solid  angle  for  scattering 
radiation  into  a  particular  direction  specified  by  the  infinitesimal 
solid  angle  dO 
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Scattered  radiation  flux  into  specified  solid  angle 
Indicent  radiation  flux 

Azimuthal  angle  of  spherical  coordinate  scattering  geometry 
shown  in  Fig.  5 

Denotes  integration  over  specified  subtended  angular  range 

Bar  over  a  quantity  denotes  average  (of  unspecified  weighting) 
oveT  the  particle  size  distribution  function  (PSDF) 

Denotes  a  set  of  parameters 


Pertains  to  quantities  with  polarization  vector  perpendicular  to 
scattering  plane 

Pertains  to  quantities  with  polarization  vector  parallel  to 
scattering  plane 

Pertains  to  directly  backscattered  radiation 
Denotes  specific  angular  scattering  direction 
Denotes  particle  chemical  species 
Pertains  to  transmitted  radiation 
Unpolarized 


Denotes  measured  quantity  after  all  calibration  corrections 
have  been  applied 

Denotes  a  particular  set  of  size  parameters,  each  element  of 
which  is  associated  with  a  given  particle  species 

Pertains  to  lower  limit  of  experimental  imprecision  band 

Pertains  to  upper  limit  of  experimental  imprecision  band 
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Denotes  a  particular  set  of  experimental  values  of  the  degree  of 
polarization,  each  element  of  which  lies  between  the  estimated 
imprecision  limits  of  the  scattering  data  from  a  given  angular 
channel 


